Solve equationsystem (A*V).'*B*(A*V)=C for matrix A
1 view (last 30 days)
Show older comments
Hello all,
im asking me now a while how can i solve this equation system:
(A*V).'*B*(A*V)=C
where V, B and C are known.
V is a vector basis. B and C are symmetric square matrices. i want to find a matrix A which is initial simply a unity square matrix. The matrix A should be like a scaling matrix for a vector basis.
Is there a possibility to solve the equation system, maybe iterative, for the matrix A?
Maybe there are more than one solution. but tell me your guesses or ideas what is maybe possible (algorithms, ideas what ever).
i have to find the minimum: min((A*V).'*B*(A*V)-C=0).
Thanks in advance and i hope for some answers :)
best regards
0 Comments
Accepted Answer
John D'Errico
on 8 Jan 2020
Edited: John D'Errico
on 8 Jan 2020
As it turns out though, your problem is actually quite simple.
Assume that V is full rank. As long as that is true, then as can arbitrarily make the substitution
X = A*V
V is known and of full rank, so if we can solve the simpler quadratic form
X'*B*X = C
then we can trivially recover A from X.
I'll note there is much to be found in the literature, solving what seems to be a similar equation, X*A*X=B. For example:
Of course, your problem is different because of the transpose. And that is what makes the solution trivial. Lets take this a step further. B is symmetric and square. Is it positive definite? If so, then we can write B as
B = L'*L
So a Cholesky decomposition of B. Now make a further transformation, with
Y = L*X
Your problem now becomes to solve for Y in the problem
Y'*Y = C
Again, is C positive definite? If so, then we can write Y in the form of a Cholesky decomposition of C. Once you have Y, recover X. Once you have X, recover A.
So the problem becomes trivial, if B and C are positive definite, and V is any full rank matrix.
Edit: Since I see by your response to Matt that V is non-square, the problem is still relatively easy, although the solution now becomes far less unique. If
X = A*V
then once X is known, recover a solution for A as
A = X*pinv(V)
7 Comments
John D'Errico
on 9 Jan 2020
It is not just that Y is not unique. The extraction of A, given X is also not unique. So there is a set of infinitely many matrices A that solve the problem.
Anyway, this is now a completely different problem. For example, suppose you had matrices B1,C1, B2,C2, with a fixed basis V? The method I showed that can generate one non-unique solution A is not directly extensible to two sets of matrices. Could you generate two solutions, A1 and A2, (or A1,A2, ... , An) then take the average? It is not clear that the average of thoe matrices, for example, would also be a solution, or even be at all meaningful.
The set of 10x2 matrix solutions for any given 2x2 matrix C, such that Y'*Y ==C is not related linearly to the matrix C. Those are essentially quadratic equations.
That leaves you with a problem where you have a 10x10 matrix of unknowns, so 100 unknowns, but still possibly not enough pieces of information to generate a unique solution.
You might decide to try a nonlinear optimizer, perhaps lsqnonlin.
More Answers (3)
Matt J
on 8 Jan 2020
Edited: Matt J
on 8 Jan 2020
for a first solution we could assume A is diagonal.
An iterative solution:
fun=@(a) objective(a,B,C,V);
a=lsqnonlin(fun, ones(length(B),1) );
A=diag(a);
function F=objective(a,B,C,V)
A=diag(a);
F=C-(A*V).'*B*(A*V);
end
8 Comments
Matt J
on 9 Jan 2020
Edited: Matt J
on 9 Jan 2020
With the modified code below, and after suitable normalization of your B,C,V data (note that this doesn't change the solutions A), I obtain a solution a that indeed is not very different from the initial guess, but it is different and measurably improves the error. As I told you, the initial guess of a=ones(N,1) was almost optimal to begin with.
First-Order Norm of
Iteration Func-count Residual optimality Lambda step
0 3331 0.000683151 0.37 0.01
1 6664 0.000681193 0.074 1 0.00211418
2 9999 0.000671141 0.0379 10000 2.57346e-05
C =
0.0042 0.0000
0.0000 227.3884
>> Error(a)
ans =
-0.0259 0.0000
0.0000 -0.0000
Modified code:
load matlab.mat
B=sparse(B/1e4);
C=C/1e14;
V=V/1e5;
Error=@(A) objective(A,B,C,V);
opts=optimoptions(@lsqnonlin,'MaxIterations',10);
a=lsqnonlin(Error, ones(length(B),1) ,[],[],opts );
function F=objective(a,B,C,V)
N=length(B);
A=spdiags(a(:),0,N,N);
F=C-(A*V).'*B*(A*V);
end
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!