Thread Subject: Solving AM = MB

I would like to solve AM = MB where A,B,M are 3x3 matrices.

What I came up with was to equate every entry in C to the corresponding entry in D (where C = AM and D = MB).

You can then vectorize M (call it Mv) and figure out the 9 equations to fill a 9x9 matrix on both sides

9x9 matrix times Mv = 9x9 matrix times Mv

Each row on the right can be subtracted from the same row on the left, leaving

9x9 times Mv = 0

The problem is, I could construct this matrix easily by hand, but this seems like an obnoxious process to write in code (a couple of loops or something?) Is there a better/different way to do this so that I can use some built in functions?

Thanks,
Dave

"David Doria" <daviddoria@gmail.com> wrote in message <gd34pk$csc$1@fred.mathworks.com>...
> I would like to solve AM = MB where A,B,M are 3x3 matrices.
>
> What I came up with was to equate every entry in C to the corresponding entry in D (where C = AM and D = MB).
>
> You can then vectorize M (call it Mv) and figure out the 9 equations to fill a 9x9 matrix on both sides
>
> 9x9 matrix times Mv = 9x9 matrix times Mv
>
> Each row on the right can be subtracted from the same row on the left, leaving
>
> 9x9 times Mv = 0
>
> The problem is, I could construct this matrix easily by hand, but this seems like an obnoxious process to write in code (a couple of loops or something?) Is there a better/different way to do this so that I can use some built in functions?


This gets asked every once in a while.
So I wrote a function to do it. There is
nice trick inside to build the necessary
matrices.

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=19614&objectType=FILE

HTH,
John

So I was looking to solve
AM = MB

And abba.m solves
AB = BA

which is slightly different, because there are only 2 matrices A and B, instead of A,B, and M.

I was a victim of your comment "% have fun with this next line...". I read the help for kron() but I have no background with tensors. Could you give a couple line explanation of what that is doing? I can't modify your code to fit this problem because I don't know what that line does, and it seems to be about the whole thing!!

Thanks,

Dave

"David Doria" <daviddoria@gmail.com> wrote in message <gd38sd$eb9$1@fred.mathworks.com>...
> So I was looking to solve
> AM = MB
>
> And abba.m solves
> AB = BA
>
> which is slightly different, because there are only 2 matrices A and B, instead of A,B, and M.
>
> I was a victim of your comment "% have fun with this next line...". I read the

Oh, heck. Victim? Thats a bit extreme. 8-)
Think of it as mathematical character building.


> help for kron() but I have no background with tensors. Could you give a couple line explanation of what that is doing? I can't modify your code to fit this problem because I don't know what that line does, and it seems to be about the whole thing!!


Yes, I'll admit that kron line is a bit crucial
to do it efficiently. It is slick though, you
gotta admit. The null call is also important.

If a solution exists, then it should take no
more than a small change. If the call to
null produces an empty matrix, then no
solution exists. Assuming that n = size(A,1),
it can be compressed into only one line.

M = reshape(mean(null(kron(eye(n),A) - ...
     kron(B',eye(n))),2),n,n);

I just have to add, "have fun with that line."
Don't ask me why, the devil on my shoulder
made me sat that. ;-)

John

"David Doria" <daviddoria@gmail.com> wrote in message <gd38sd$eb9$1@fred.mathworks.com>...
> So I was looking to solve
> AM = MB
>
> And abba.m solves
> AB = BA
>
> which is slightly different, because there are only 2 matrices A and B, instead of A,B, and M.
>
> I was a victim of your comment "% have fun with this next line...". I read the help for kron() but I have no background with tensors. Could you give a couple line explanation of what that is doing? I can't modify your code to fit this problem because I don't know what that line does, and it seems to be about the whole thing!!
>
> Thanks,
>
> Dave


John's method is difficult to generalize for all cases. I have tried to generalize his technique for my own uses.

Assuming you have a matrix equation of the form:

  A1*X + X*A2 + A3*X*A4 = 0

where A1, A2, A3, A4, and X are all n-by-n matrices, you can solve this equation using the following commands:

A1X = kron(eye(n),A1);
XA2 = kron(A2.',eye(n));
A3XA4 = kron(eye(n),A3)*kron(A4.',eye(n));
M = (A1X + XA2 + A3XA4);
Mnull = null(M);
if ~isempty(Mnull)
  X = reshape(sum(Mnull,2),n,n);
else
  X = zeros(n,n);
end

However, if you have a matrix equation of the form:

A1*X + X*A2 + A3*X*A4 = A5

I use the variation:

% Same here
A1X = kron(eye(n),A1);
XA2 = kron(A2.',eye(n));
A3XA4 = kron(eye(n),A3)*kron(A4.',eye(n));
M = (A1X + XA2 + A3XA4);
% Here's where it differs
x = M\A5(:);
X = reshape(x,n,n);

There are improvements which can be made here, but I hope this gives you the gist of it.

I'm still have no idea what that is doing!

I'm not (yet) concerned with that general case of multiple combinations of the matrix multiplied by other matrices. What I would like to know though is if my original explanation of what is going on is the same thing that you guys are doing with this big tensor kron() thing?

Also - can someone point me to a simple, practical, engineering approach guide to tensors (ie. something that may help me understand this bit of code)? I've always heard them thrown around as a word people say when they start to sound like they don't really know what they're talking about! haha clearly you guys do not fit that bill, so you had to have learned it somewhere, maybe you can share?

Thanks!
Dave

"Tim Farajian" <timfarajian_REMOVEME@socal.rr.com.ANDREMOVEME> wrote in message <gd3i4i$13v$1@fred.mathworks.com>...
> John's method is difficult to generalize for all cases. I have tried to generalize his technique for my own uses.

(snip)

> There are improvements which can be made here, but I hope this gives you the gist of it.

Nicely done Tim.

John

"David Doria" <daviddoria@gmail.com> wrote in message <gd3jla$d35$1@fred.mathworks.com>...
> I'm still have no idea what that is doing!
>
> I'm not (yet) concerned with that general case of multiple combinations of the matrix multiplied by other matrices. What I would like to know though is if my original explanation of what is going on is the same thing that you guys are doing with this big tensor kron() thing?
>
> Also - can someone point me to a simple, practical, engineering approach guide to tensors (ie. something that may help me understand this bit of code)? I've always heard them thrown around as a word people say when they start to sound like they don't really know what they're talking about! haha clearly you guys do not fit that bill, so you had to have learned it somewhere, maybe you can share?


Dave,

Ok, I'll show some mercy. It is rare, but every
once in a while...

Wikipedia has a reasonably good description,
and it is convenient to find.

http://en.wikipedia.org/wiki/Kronecker_product

Down towards the bottom, under matrix equations,
you will find the claim that the matrix product
A*X*B can be written in terms of the unrolled
form of X as a vector. In MATLAB form, it reduces
to:

  kron(A,B.')*X(:)

You can also express the simpler forms A*X or
X*B using the above, by putting in an nxn identity
matrix for one or the other elements in the
expression.

Given that, everything else should now fall into
place.

John

"David Doria" <daviddoria@gmail.com> wrote in message <gd3jla$d35$1@fred.mathworks.com>...
> I'm still have no idea what that is doing!
>
> I'm not (yet) concerned with that general case of multiple combinations of the matrix multiplied by other matrices. What I would like to know though is if my original explanation of what is going on is the same thing that you guys are doing with this big tensor kron() thing?
>
> Also - can someone point me to a simple, practical, engineering approach guide to tensors (ie. something that may help me understand this bit of code)? I've always heard them thrown around as a word people say when they start to sound like they don't really know what they're talking about! haha clearly you guys do not fit that bill, so you had to have learned it somewhere, maybe you can share?
>
> Thanks!
> Dave

Maybe a short explanation to get you started:

John used the word "unrolling". By this, he means he considers the (for example) 3-by-3 matrix as a 9-by-1 vector. So the matrix X:

[x11 x12 x13]
[x21 x22 x23]
[x31 x32 x33]

becomes the vector x:

[x11]
[x21]
[x31]
[x12]
[x22]
[x32]
[x13]
[x23]
[x33]

Then, the operation A*X is equivalent to then operation A2*x where

A2 = kron(eye(n),A);

Explore the form of A2 and see if you can figure out why.

The vector resulting from A2*x then must be reshaped back into a 3-by-3 matrix to obtain A*X.

For right multiplication, X*B is the equivalent to B2*x where

B2 = kron(B.',eye(n));

Again, explore the form of B2 to see if you understand why.

Now,
 
  A1*X-X*A2 = 0 % (3-by-3) matrices

becomes

  A2*x-B2*x = (A2-B2)*x = 0 % (9-by-1) vectors

All you must do here is solve for the vector x and "reshape" it back into a 3-by-3 matrix.


> so you had to have learned it somewhere, maybe you can
> share?


Not necessarily true. I learned this method from John, but John is one of those men who doesn't necessarily learn it from anywhere. He invents.

kron is neat. Here is another way:

% Data
A=eye(3); A(3)=1;
B=eye(3); B(3)=1;

n=length(A);

% Engine
c(1:n)={A};
BIGA = blkdiag(c{:});

c(1:n)={B.'};
BIGB = blkdiag(c{:});
k=reshape(1:n*n,n,n)';
BIGB=BIGB(k(:),k(:));

M = null(BIGA-BIGB);
% All combination of M(n,n,:) is solution
k=size(M,2);
M = reshape(M,[n n k])

% Test
% random combination
v=rand(1,1,k);
Mrand=bsxfun(@times,M,v);
Mrand=sum(Mrand,3)
% Check
A*Mrand
Mrand*B

% Bruno

Hi all,

This program is interesting. But would the solution be unique?? I have tried the algorithms presented here, but it does not give me a solution that makes sense (because I have the known solution). However, checking XA and BX, they are indeed equal. So it makes me think that the solution is not unique. Any inputs here?

"Pinpress" <nothing@nothing.edu> wrote in message <gg9pp2$q6m$1@fred.mathworks.com>...

> However, checking XA and BX, they are indeed equal. So it makes me think that the solution is not unique. Any inputs here?

That why I wrote:
% All combination of M(n,n,:) is solution

Bruno

OK. Let me provide an example:

A =

    0.9932 -0.1156 0.0144 27.0156
    0.1150 0.9926 0.0385 0.5268
   -0.0187 -0.0366 0.9992 4.8570
         0 0 0 1.0000
B =

    0.9932 0.1160 0.0141 -41.1405
   -0.1154 0.9926 -0.0374 5.7466
   -0.0183 0.0355 0.9992 -8.3134
         0 0 0 1.0000

I have followed the algorithm to compute X: X*A = B*X, and obtained:

X =

   -0.5700 0.0488 -0.0716 0.9674
    0.0654 0.5853 -0.0091 -0.0883
   -0.0465 0.0426 -0.4048 1.4363
    0.0000 -0.0000 0.0000 0.3822

But my "real" solution, which is known, is:

T =
   -0.9999 0.0098 -0.0106 44.6588
    0.0098 1.0000 0.0001 127.7536
    0.0106 -0.0000 -0.9999 -30.5660
         0 0 0 1.0000

You can verify that X*A = B*X, as well as T*A = B*T. By the way, I really need a orthogonal matrix as my solution. So, X is NOT what I want.

>> det(X)

ans =

    0.0512

>> det(T)

ans =

    1.0000


"Pinpress" <nothing@nothing.edu> wrote in message <gg9pp2$q6m$1@fred.mathworks.com>...
> Hi all,
>
> This program is interesting. But would the solution be unique?? I have tried the algorithms presented here, but it does not give me a solution that makes sense (because I have the known solution). However, checking XA and BX, they are indeed equal. So it makes me think that the solution is not unique. Any inputs here?

I see -- sorry I overlooked.

So here is my real problem to solve:

I have many equations like:

X*A1 = B1*X;
X*A2 = B2*X;
 etc.

And I know X is orthogonal, because it is a rigid transformation matrix. How do I construct or recover a sensible X from these equations? Thanks for any input.

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gg9ql5$6va$1@fred.mathworks.com>...
> "Pinpress" <nothing@nothing.edu> wrote in message <gg9pp2$q6m$1@fred.mathworks.com>...
>
> > However, checking XA and BX, they are indeed equal. So it makes me think that the solution is not unique. Any inputs here?
>
> That why I wrote:
> % All combination of M(n,n,:) is solution
>
> Bruno

"Pinpress" <nothing@nothing.edu> wrote in message <gg9qui$9si$1@fred.mathworks.com>...
> I see -- sorry I overlooked.
>
> So here is my real problem to solve:
>
> I have many equations like:
>
> X*A1 = B1*X;
> X*A2 = B2*X;
> etc.
>
> And I know X is orthogonal, because it is a rigid transformation matrix. How do I construct or recover a sensible X from these equations? Thanks for any input.
>

You can stack many of these on top of each other, solving
for one overall matrix.

But you cannot build in the requirement of orthogonality into
a linear least squares solution.

John

Sorry John, I did not quite understand.

All my matrices, X, Ai, Bi, are 4-by-4 square matrix. My gut feeling is that, although for each X*A = B*X equation, there many possible X's to satisfy, but there should be a unique one that satisfies all the equations because all the equations are derived from physical quantities. Then I am kind stuck on how to find it.

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gg9sua$4gv$1@fred.mathworks.com>...
> "Pinpress" <nothing@nothing.edu> wrote in message <gg9qui$9si$1@fred.mathworks.com>...
> > I see -- sorry I overlooked.
> >
> > So here is my real problem to solve:
> >
> > I have many equations like:
> >
> > X*A1 = B1*X;
> > X*A2 = B2*X;
> > etc.
> >
> > And I know X is orthogonal, because it is a rigid transformation matrix. How do I construct or recover a sensible X from these equations? Thanks for any input.
> >
>
> You can stack many of these on top of each other, solving
> for one overall matrix.
>
> But you cannot build in the requirement of orthogonality into
> a linear least squares solution.
>
> John

"Pinpress" <nothing@nothing.edu> wrote in message <gg9th2$bo1$1@fred.mathworks.com>...
> Sorry John, I did not quite understand.
>
> All my matrices, X, Ai, Bi, are 4-by-4 square matrix. My gut feeling is that, although for each X*A = B*X equation, there many possible X's to satisfy, but there should be a unique one that satisfies all the equations because all the equations are derived from physical quantities. Then I am kind stuck on how to find it.
>

Look back at this thread. We showed how to solve
ONE problem, of the required form, for a 4x4 matrix.
This is transformed into a system of equations in 16
unknowns.

Now, simply concatenate the corresponding matrices
vertically on top of each other, solving one large
problem. This will give a global solution to all problems,
if one exists. It will still not be unique, nor need it be
an orthogonal matrix.

John