How to find matrix S, given by this equation
5 views (last 30 days)
John D'Errico on 24 Apr 2022
There is no unique solution. Or, there may be no solution at all.
This looks vaguely like an algebraic Riccati equation. Your problem is of the general form
inv(S)*L*S = J
Where S is the unknown matrix. If the matrix S has an inverse, then L must have the same rank as your right hand side matrix. And we know that the right hand side matrix (J) is rank deficient.
First, make up a random problem.
J = zeros(4,4);
Jdel = magic(3);
J(2:end,2:end) = Jdel
So the right hand side is a matrix of the indicated form. Now I'll choose a known non-singular matrix S. Here is a simple one:
S0 = toeplitz(1:4)
L = S0*J*inv(S0)
So the above problem must satisfy the equation: inv(S)*L*S==J. A solution MUST exist.
Now the problem becomes, given only J and L, can we recover S? I'll pretend I do not know the matrix S0.
Can we find a solution of the form L*S = S*J? Note that I have implicitly left multiplied by S to arrive at that form. But as long as S is non-singular, that is legal, and since you want a solution where S has an inverse, then S MUST be non-singular.
The solution is not too nasty. It uses what I long ago called the kron trick.
n = size(J,1);
M = kron(eye(n),L) - kron(J.',eye(n));
Note that M will be a singular matrix. The problem is, in this case M will have rank n^2-n, not just n^2-1.
We need to solve a homogeneous linear system now.
Mnull = null(M)
ANY linear combination of the columns of Mnull will suffice. So I will choose some random linear combination.
S = reshape(Mnull*randn(n,1),n,n)
Does this solve your problem? Well, yes. Sort of. But it is not a unique solution to that problem in any form.
norm(inv(S)*L*S - J)
Which results in floating point trash. So we have recovered a solution, but not the one I started with. The point is, again, there is no unique solution. And had I not been careful in how I created L then no solution at all may have existed.
More Answers (2)
Torsten on 24 Apr 2022
L is NxN and 0_(N-1) is (N-1)x1 ? J_delta is a full (N-1)x(N-1) matrix ?
Then use the eigenvector corresponding to the eigenvalue 0 of L to build S.
Sam Chak on 24 Apr 2022
@NA, thanks for your clarification on the dimensions.
Is this some kind of a Linear Algebra problem?
It seems that the square matrix L is unique, and has to satisfy certain conditions.
For example, say the matrix is given by
L = [0 1 0; 0 0 1; 0 -6 -5]
0 1 0
0 0 1
0 -6 -5
Then, the eigenvalues of matrix are
s = eig(L)
Next, the matrix is constructed from the eigenvalues such that
S = [1 1 1; s(1) s(2) s(3); s(1)^2 s(2)^2 s(3)^2]
1.0000 1.0000 1.0000
0 -2.0000 -3.0000
0 4.0000 9.0000
Finally, we can test if is a diagonal matrix where the diagonal elements are the eigenvalues of matrix :
D = S\L*S
0 0.0000 0.0000
0 -2.0000 -0.0000
0 0.0000 -3.0000
Only a centain square matrix can produce such diagonal matrix. For example:
L = [0 1 0; 0 0 1; 0 -15 -8]
0 1 0
0 0 1
0 -15 -8