MATLAB Central Newsreader - finding a matrix used in matrix transformation

finding a matrix used in matrix transformation

pushkarini — Tue, 09 Nov 2010 17:46:05 +0000

Hi
I have known matrices A and B (both square of size n say) such that
B = inv(M)*A*M
so that M*B = A*M
I need to find M
I can solve a system of n^2 equations in n^2 unknowns to find M, but is there an easier matlab command or another trick to do this?

thanks

Re: finding a matrix used in matrix transformation

James Tursa — Tue, 09 Nov 2010 17:58:04 +0000

Re: finding a matrix used in matrix transformation

Matt J — Tue, 09 Nov 2010 18:06:04 +0000

"pushkarini " <pushkarini.a@gmail.com> wrote in message <ibc1cs$61o$1@fred.mathworks.com>...
> Hi
> I have known matrices A and B (both square of size n say) such that
> B = inv(M)*A*M
> so that M*B = A*M
=====

Do you know for a fact that a nonsingular solution M exists?

Let m=M(:). and
Q=(kron(B.', speye(n)) - kron(speye(n),A))

Then the equation M*B-A*M is equivalent to

Q*m=0

So, if Q is nonsingular, the only solution is m=0, corresponding to M=zeros(n)
which is singular. This is why I asked if you know for a fact that a nonsingular solution M exists.

If you're sure that Q is singular, then you should be able to recover m as an eigenvector of Q with eigenvalue 0. See "doc eig"

Re: finding a matrix used in matrix transformation

Bruno Luong — Tue, 09 Nov 2010 18:42:03 +0000

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ibc2ic$mjm$1@fred.mathworks.com>...

> If you're sure that Q is singular, then you should be able to recover m as an eigenvector of Q with eigenvalue 0. See "doc eig"

Or simply using NULL which gives a basis of the null space (at the accuracy of the floating point).

Bruno

Re: finding a matrix used in matrix transformation

Roger Stafford — Tue, 09 Nov 2010 19:40:04 +0000

"pushkarini " <pushkarini.a@gmail.com> wrote in message <ibc1cs$61o$1@fred.mathworks.com>...
> Hi
> I have known matrices A and B (both square of size n say) such that
> B = inv(M)*A*M
> so that M*B = A*M
> I need to find M
> I can solve a system of n^2 equations in n^2 unknowns to find M, but is there an easier matlab command or another trick to do this?
>
> thanks
- - - - - - - - - -
  The requirement that an invertible matrix M exists such that B = inv(M)*A*M is an equivalence relation and A and B are called 'similar' matrices if it holds between them. It has been shown that matrices are similar if and only if their Jordan normal forms are the same up to possible rearrangement of the Jordan blocks. There is a good paper on this subject by Kahan at:

http://www.eecs.berkeley.edu/~wkahan/MathH110/jordan.pdf

  Presumably the way to find the M you seek would be to convert by similarity transform both A and B to their respective Jordan normal forms. If these are alike up to rearrangement of their Jordan blocks, then combining all of these transformations and permutations would give you the necessary M matrix. If not, then no such M can exist.

  Unfortunately I don't know enough about the subject to tell you how to find the Jordan normal form using matlab, but reading the above paper might give you some ideas. From a numerical analysis point of view finding the Jordan normal form is known to be numerically unstable in that small errors can make large differences in the result. However I suspect this probably reflects the inherent difficulty of solving the problem you pose.

Roger Stafford

Re: finding a matrix used in matrix transformation

Roger Stafford — Tue, 09 Nov 2010 22:04:04 +0000

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ibc82k$ui$1@fred.mathworks.com>...
> "pushkarini " <pushkarini.a@gmail.com> wrote in message <ibc1cs$61o$1@fred.mathworks.com>...
> > Hi
> > I have known matrices A and B (both square of size n say) such that
> > B = inv(M)*A*M
> > so that M*B = A*M
> > I need to find M
> > I can solve a system of n^2 equations in n^2 unknowns to find M, but is there an easier matlab command or another trick to do this?
> >
> > thanks
> - - - - - - - - - -
> The requirement that an invertible matrix M exists such that B = inv(M)*A*M is an equivalence relation and A and B are called 'similar' matrices if it holds between them. It has been shown that matrices are similar if and only if their Jordan normal forms are the same up to possible rearrangement of the Jordan blocks. There is a good paper on this subject by Kahan at:
>
> http://www.eecs.berkeley.edu/~wkahan/MathH110/jordan.pdf
>
> Presumably the way to find the M you seek would be to convert by similarity transform both A and B to their respective Jordan normal forms. If these are alike up to rearrangement of their Jordan blocks, then combining all of these transformations and permutations would give you the necessary M matrix. If not, then no such M can exist.
>
> Unfortunately I don't know enough about the subject to tell you how to find the Jordan normal form using matlab, but reading the above paper might give you some ideas. From a numerical analysis point of view finding the Jordan normal form is known to be numerically unstable in that small errors can make large differences in the result. However I suspect this probably reflects the inherent difficulty of solving the problem you pose.
>
> Roger Stafford
- - - - - - - - - -
I just found out that the Symbolic Toolbox has a function, 'jordan', which will find the Jordan normal form, either symbolically or numerically, of any square matrix and the corresponding similarity transform. Therefore this, together with my previous remarks, should solve your problem.

Roger Stafford

Re: finding a matrix used in matrix transformation

Matt J — Tue, 09 Nov 2010 23:02:04 +0000

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ibc4lr$cjv$1@fred.mathworks.com>...
> "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ibc2ic$mjm$1@fred.mathworks.com>...
>
> > If you're sure that Q is singular, then you should be able to recover m as an eigenvector of Q with eigenvalue 0. See "doc eig"
>
> Or simply using NULL which gives a basis of the null space (at the accuracy of the floating point).
========

This seems to work, as tested below. Unfortunately, though, null() won't work for sparse Q, so it tends to be quite slow to apply it to full Q

n=15;

%simulated data
Mtrue=rand(n);
A=rand(n);
B=inv(Mtrue)*A*Mtrue;

%Engine
Q=full( kron(B.', speye(n)) - kron(speye(n),A) );

N=null(Q);

d=size(N,2);
dets=zeros(1,d);
for ii=1:d

   M=reshape(N(:,ii),n,n);

  dets(ii)=abs(det(M));

end

[~,j]=max(dets);

M=reshape(N(:,ii),n,n);

Re: finding a matrix used in matrix transformation

Bruno Luong — Tue, 09 Nov 2010 23:20:04 +0000

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ibcjtb$1nt$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ibc4lr$cjv$1@fred.mathworks.com>...
> > "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ibc2ic$mjm$1@fred.mathworks.com>...
> >
> > > If you're sure that Q is singular, then you should be able to recover m as an eigenvector of Q with eigenvalue 0. See "doc eig"
> >
> > Or simply using NULL which gives a basis of the null space (at the accuracy of the floating point).
> ========
>
> This seems to work, as tested below. Unfortunately, though, null() won't work for sparse Q, so it tends to be quite slow to apply it to full Q

That's why I have submitted a null for sparse matrix in FEX. That requires recent Matlab version thought.

I'm still thinking about Roger's Jordan form which should be provides a solution of the eqt

B = inv(M)*A*M

which is not obvious to get from

the solutions of

M*B = A*M, (provides by the kron/null for example)

I can't find a simple example where Jordan method gives a more focus solution space. Does someone can build of such example?

Bruno

Re: finding a matrix used in matrix transformation

Roger Stafford — Wed, 10 Nov 2010 03:30:10 +0000

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ibckv4$8q2$1@fred.mathworks.com>...
> ........
> I'm still thinking about Roger's Jordan form which should be provides a solution of the eqt
>
> B = inv(M)*A*M
>
> which is not obvious to get from
>
> the solutions of
>
> M*B = A*M, (provides by the kron/null for example)
>
> I can't find a simple example where Jordan method gives a more focus solution space. Does someone can build of such example?
>
> Bruno
- - - - - - - -
  I would see it proceeding as follows:

[V1,J1] = jordan(A);
[V2,J2] = jordan(B);

If J1 and J2 are alike except for a rearrangement of Jordan blocks, let P be the permutation matrix that makes J1 and J2 match blocks, so that J2 = inv(P)*J1*P. Then setting

M = V1*P*inv(V2);

gives the required M. That is because

inv(M)*A*M = (V2*inv(P)*inv(V1)) * A * (V1*P*inv(V2)) =
V2*inv(P)*J1*P*inv(V2) = V2*J2*inv(V2) = B.

  The step that takes additional work is finding the necessary permutation that makes J1 and J2 match.

  I do worry about the sensitivity of 'jordan' to small rounding errors. It could conceivably cause J1 and J2 to have different Jordan blocks when they ought to be alike.

  To give the "kron/null" method and this one a thorough test one could construct a Jordan normal form matrix with a set of a number of Jordan blocks which all have, say, equal eigenvalues. (Perhaps these could even be all zero.) Then with two random but invertible matrices do similarity transforms on the normal matrix to produce A and B. Then try out the kron/null and the above methods to see if either has troubles.

Roger Stafford

Re: finding a matrix used in matrix transformation

Bruno Luong — Wed, 10 Nov 2010 07:27:07 +0000

Roger,

Some preliminary facts:

Let
problem 1, B = inv(M)*A*M
problem 2, M*B = A*M

Any solution of problem 1 is solution of problem 2. So space-solution of problem 2 contains space solution of problem 1.

I don't have symbolic toolbox so I can't check Jordan method. However if I generate A and B such as required, then I check the Kron/Null and it gives not 1-dimensional space solution but n-dimensional space-solution (this should be provable by Kronecker property, but right now I just observe this fact).

What I wonder is this:

In the problem 1, we look for the Jordan-transformation matrix M (from A to B) that is not singular. But could the KRON/NULL basis contains somehow singular matrices, yet by combining them we get non-singular solution? If yes I would like to see an example.

Here is the code I use to check Kron/Null

% Generate matrices
J1 = [1 1 0; 0 1 0; 0 0 2];
I = eye(3);
P = I(:,[2 3 1]);
J2 = inv(P)*J1*P;
V1 = ceil(99*rand(3,3));
V2 = ceil(99*rand(3,3));

A = V1*J1*inv(V1)
B = V2*J2*inv(V2)
% Jordan transformation
M = V1*P*inv(V2);

fprintf('Jordan check\n');
fprintf(' M*B:\n');
disp(M*B)
fprintf(' A*M:\n');
disp(A*M)

%% Kron/Null
I = eye(size(A));
Q = kron(B.',I)-kron(I,A);
N = null(Q);
% Try to recover M from N
x = N\M(:);
Mrecov = reshape(N*x,size(A));
errorn = norm(Mrecov-M)/norm(M);
if errorn < 1e-10
    fprintf('OK: We can recover Jordan solution\n')
end
% Check if each basis matrix is the solution
for k=1:size(N,2)
    fprintf('Kron/NULL solution #%d check\n', k);
    M = reshape(N(:,1),[3 3]);
    if rank(M) < length(A)
        fprintf('Kron/NULL method only solve the M*B=A*M\n');
    else
        fprintf('Kron/NULL method solve the harder problem B=inv(M)*A*M\n');
    end
    fprintf(' M*B:\n');
    disp(M*B)
    fprintf(' A*M:\n');
    disp(A*M)
end

% Bruno

Re: finding a matrix used in matrix transformation

pushkarini — Wed, 10 Nov 2010 08:31:03 +0000

Thanks a lot fo your quick reply!

The method using kron command works. the one using jordan matrices gives me the error
Error: Similarity matrix too large [mljordan], also it took a loong time to solve

Re: finding a matrix used in matrix transformation

Matt J — Wed, 10 Nov 2010 15:20:05 +0000

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ibdhgb$k7j$1@fred.mathworks.com>...
>
> What I wonder is this:
>
> In the problem 1, we look for the Jordan-transformation matrix M (from A to B) that is not singular. But could the KRON/NULL basis contains somehow singular matrices, yet by combining them we get non-singular solution? If yes I would like to see an example.
======

You mean if the KRON/NULL basis contains *only* singular matrices? That would indeed be a problem.

To verify this, I guess you would want to solve a max-min squares problem. In particular, if M(:,:,i) are KRON/NULL basis matrices, then you want to find
coefficients c(i) and
an n x 1 unit vector u satisfying

  norm( sum_i c(i)*M(:,:,i)*u)>0

which you might do by solving

  max( min ( norm( sum_i c(i)*M(:,:,i)*u)^2 ))

  s.t. norm(u)^2=1

where the max is over the c(i) and the min is over u.

This stretches my optimization know-how, but I guess you could might be able to do it iteratively by alternatingly maximizing over c and minimizing over u, both of which are simple computations....

Re: finding a matrix used in matrix transformation

Matt J — Wed, 10 Nov 2010 15:34:04 +0000

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ibed75$jiu$1@fred.mathworks.com>...
>
> To verify this, I guess you would want to solve a max-min squares problem. In particular, if M(:,:,i) are KRON/NULL basis matrices, then you want to find
> coefficients c(i) and
> an n x 1 unit vector u satisfying
>
> norm( sum_i c(i)*M(:,:,i)*u)>0
=====

Sorry, this should have read, find c(i) such that

min_u norm(sum_i c(i)*M(:,:,i)*u)/norm(u) >0

Re: finding a matrix used in matrix transformation

Bruno Luong — Wed, 10 Nov 2010 15:46:04 +0000

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message
> You mean if the KRON/NULL basis contains *only* singular matrices? That would indeed be a problem.

Matt,

My question is not really how to derive a non-singular M from the null space vectors N. I'm not there yet!!

My question is simpler: can it happens? Since this is what could be a drawback where Jordan approach does not suffer. To me it is not clear the motivation of Roger recommends for Jordan's transformation.

Bruno

Re: finding a matrix used in matrix transformation

Matt J — Wed, 10 Nov 2010 15:53:04 +0000

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ibdhgb$k7j$1@fred.mathworks.com>...
>
>But could the KRON/NULL basis contains somehow singular matrices, yet by combining them we get non-singular solution? If yes I would like to see an example.
======

The simplest example is probably when A=B=Mtrue=eye(n).

In that case the Kron-derived matrix is Q=zeros(n^2)
and null(Q)=eye(n^2).

In this case, each and every null matrix M derived from reshaping a column of null(Q)
is singular.

Observe also that, in addition to M=eye(n), every possible nxn nonsingular matrix M is in the span of null(Q), and is a valid solution to problem 1.

Re: finding a matrix used in matrix transformation

Bruno Luong — Wed, 10 Nov 2010 16:01:06 +0000

> The simplest example is probably when A=B=Mtrue=eye(n).
>
> In that case the Kron-derived matrix is Q=zeros(n^2)
> and null(Q)=eye(n^2).
>
> In this case, each and every null matrix M derived from reshaping a column of null(Q)
> is singular.
>
> Observe also that, in addition to M=eye(n), every possible nxn nonsingular matrix M is in the span of null(Q), and is a valid solution to problem 1.

OK. That works. Now the question is what Jordan method gives in this case?

Bruno

Re: finding a matrix used in matrix transformation

pushkarini — Wed, 10 Nov 2010 16:36:03 +0000

Hi
The kron command works for a matrix A and another one B (which i derived from A using some column and row operations.)
Next, I used orth(A) command to get an orthonormal basis for A. I then did:
B = orth(A)
and tried to find an M such that M*B = A*M, and after using the method involving kron got some value for M. But, when I used this M to compute my B back (using B = inv(M)*A*M), I did not get my B back. Why is this so? Then, I thought, it is probably not possible to convert A into orth(A) using such a method in the first place? I guess, adding proper 'lengths' to the row vectors of the orthonormal basis will solve the problem?
I am doing all this to convert my matrix A which I start with to a normal matrix (ie. whose eigenvectors are perpendicular to eachother).

Re: finding a matrix used in matrix transformation

Matt J — Wed, 10 Nov 2010 17:08:05 +0000

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ibefk2$2hn$1@fred.mathworks.com>...
>
> OK. That works. Now the question is what Jordan method gives in this case?
=====

It will give you some non-singular matrix M. Which one depends on which of the infinite possible Jordan decompositions of A=eye(n) and B=eye(n) you decide to use.

Re: finding a matrix used in matrix transformation

Matt J — Wed, 10 Nov 2010 17:36:05 +0000

"pushkarini " <pushkarini.a@gmail.com> wrote in message <ibehlj$ijd$1@fred.mathworks.com>...

Then, I thought, it is probably not possible to convert A into orth(A) using such a method in the first place?
==============

Correct. The transformation B=inv(M)*A*M will cause B to have the same eigenvalues as A. But orth(A) in general will not have the same eigenvalues as A.

I guess, adding proper 'lengths' to the row vectors of the orthonormal basis will solve the problem?
============

I don't know what this means, but it sounds doubtful.

> I am doing all this to convert my matrix A which I start with to a normal matrix (ie. whose eigenvectors are perpendicular to eachother).
====================

Always a good idea to say what you really want in the first place!!!

Assuming A has a full set of n eigenvalues, you can do the following

%Engine
[V,D]=eig(A);
iV=inv(V);

M=(V*V')^(.5);

B=inv(M)*A*M;

Re: finding a matrix used in matrix transformation

Bruno Luong — Wed, 10 Nov 2010 18:22:03 +0000

>
> max( min ( norm( sum_i c(i)*M(:,:,i)*u)^2 ))
>

Note that the above criteria (minimum singular value) is not helpful since max -> Inf if c=s*c' for s-> inf given c' any arbitrary coefficients such that sum_i c'(i)*M(:,:,i) is non singular.

It is better to minimize some quantity related to condition number or stretch of the spectrum.

But anyway the problem of optimizing eigen/singular values is slightly OT.

Bruno

Re: finding a matrix used in matrix transformation

Matt J — Wed, 10 Nov 2010 18:42:04 +0000

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ibensb$8r6$1@fred.mathworks.com>...
>
> >
> > max( min ( norm( sum_i c(i)*M(:,:,i)*u)^2 ))
> >
>
> Note that the above criteria (minimum singular value) is not helpful since max -> Inf if c=s*c' for s-> inf given c' any arbitrary coefficients such that sum_i c'(i)*M(:,:,i) is non singular.
======

I proably meant to have norm(c) bounded to 1 as well as norm(u).

In any case, if you're going to solve this iteratively, it won't matterif max -> Inf. As soon as you reach an iteration where

M=sum_i c(i)*Mbasis(:,:,i)

is non-singular, you would terminate the algorithm, with M as your solution.

Re: finding a matrix used in matrix transformation

Bruno Luong — Wed, 10 Nov 2010 19:39:04 +0000

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <ibep1s$pob$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ibensb$8r6$1@fred.mathworks.com>...
> >
> > >
> > > max( min ( norm( sum_i c(i)*M(:,:,i)*u)^2 ))
> > >
> >
> > Note that the above criteria (minimum singular value) is not helpful since max -> Inf if c=s*c' for s-> inf given c' any arbitrary coefficients such that sum_i c'(i)*M(:,:,i) is non singular.
> ======
>
> I proably meant to have norm(c) bounded to 1 as well as norm(u).
>

Yes. norm(c)=1 is equivalent to Frobenious norm of transition matrix = 1, which is bounded by the spectral norm (both norms are equivalent). So by maximizing the smallest singular-value with |c|=1, it's close to the problem of minimizing the condition number without imposing constraint on c.

> In any case, if you're going to solve this iteratively, it won't matterif max -> Inf. As soon as you reach an iteration where
>
> M=sum_i c(i)*Mbasis(:,:,i)
>
> is non-singular, you would terminate the algorithm, with M as your solution.

Nope, that's not satisfying justification IMO.

Bruno