5x5 matrix (P) with all 25 values unknown. I have two known matrices of size 5x5 each (A1 and A2) and the relation P*A1*inv(P)=A2. How do I obtain P?

A2=[1.1581 0.0279 -0.0739 -0.0169 -0.0100

-0.9172 0.7637 0.5688 0.1299 0.0795

12.8599 2.5037 -5.4500 -1.4719 -0.8824

-39.7814 -7.8095 20.0738 5.5809 2.7482

-27.9052 -5.5369 14.1918 3.2387 2.9446];

A1=[ -1.0093 0.0645 -0.0795 -0.0172 -0.5886

4.1505 0.0040 -0.6850 -0.0824 -2.9194

1.1250 0.7191 -2.4753 -0.0098 -0.3402

-1.1250 1.1521 -0.2040 -1.2556 -0.4159

2.5296 0.0455 -0.1234 -0.0646 -3.7925];

syms P a1 b1 c1 d1 e1 a2 b2 c2 d2 e2 a3 b3 c3 d3 e3 a4 b4 c4 d4 e4 a5 b5 c5 d5 e5

P1=[a1 b1 c1 d1 e1

a2 b2 c2 d2 e2

a3 b3 c3 d3 e3

a4 b4 c4 d4 e4

a5 b5 c5 d5 e5];

V=vpa(P1*A1*inv(P1),6);

solve(V==A2)

%%This gives all 0s.

%%I also tried the below code but still getting zeroes.

global A1 A2 B1 B2

guess=[10;0.1; 0.1; 0.1; 0.1; 0.1; 10; 0.1; 0.1; 0.1; 0.1; 0.1; 10; 0.1; 0.1; 0.1; 0.1; 0.1; 1; 0.1; 0.1; 0.1; 0.1; 0.1; 1];

%guess=[1 -1 -2 1 2 1 1 -3 2 4 -2 1 1 3 2 3 -3 1 1.5 2.5 2.8 3 3.5 3.2 1];

%guess=[1 0.1 0.1 0.1 0.1 0.1 1 0.1 0.1 0.1 0.1 0.1 1 0.1 0.1 0.1 0.1 0.1 1 0.1 0.1 0.1 0.1 0.1 1];

options = optimoptions(@lsqnonlin,'OptimalityTolerance', 1e-16, 'FunctionTolerance', 1e-16,'StepTolerance',1e-16);

result=lsqnonlin(@eqns,guess,[],[],options);

%result=fsolve(@eqns,guess,options);

P=[result(1) result(6) result(11) result(16) result(21)

result(2) result(7) result(12) result(17) result(22)

result(3) result(8) result(13) result(18) result(23)

result(4) result(9) result(14) result(19) result(24)

result(5) result(10) result(15) result(20) result(25)];

function fncs = eqns(z)

global A1 A2 B1 B2

a1=z(1);

a2=z(2);

a3=z(3);

a4=z(4);

a5=z(5);

b1=z(6);

b2=z(7);

b3=z(8);

b4=z(9);

b5=z(10);

c1=z(11);

c2=z(12);

c3=z(13);

c4=z(14);

c5=z(15);

d1=z(16);

d2=z(17);

d3=z(18);

d4=z(19);

d5=z(20);

e1=z(21);

e2=z(22);

e3=z(23);

e4=z(24);

e5=z(25);

P1=[a1 b1 c1 d1 e1

a2 b2 c2 d2 e2

a3 b3 c3 d3 e3

a4 b4 c4 d4 e4

a5 b5 c5 d5 e5];

expres=P1*A1*inv(P1)-A2;

fncs(1)= expres(1,1);

fncs(2)= expres(1,2);

fncs(3)= expres(1,3);

fncs(4)= expres(1,4);

fncs(5)= expres(1,5);

fncs(6)= expres(2,1);

fncs(7)= expres(2,2);

fncs(8)= expres(2,3);

fncs(9)= expres(2,4);

fncs(10)= expres(2,5);

fncs(11)= expres(3,1);

fncs(12)= expres(3,2);

fncs(13)= expres(3,3);

fncs(14)= expres(3,4);

fncs(15)= expres(3,5);

fncs(16)= expres(4,1);

fncs(17)= expres(4,2);

fncs(18)= expres(4,3);

fncs(19)= expres(4,4);

fncs(20)= expres(4,5);

fncs(21)= expres(5,1);

fncs(22)= expres(5,2);

fncs(23)= expres(5,3);

fncs(24)= expres(5,4);

fncs(25)= expres(5,5);

fncs(26)= expres(4,5);

end

Even if we assume that A1 and A2 do not satisfy this relation A2=P*A1*inv(P), how do we find P for some other A1 A2 which do satisfy this relation. Please help!

[Richa Dubey]

Lets take A1 and A2 both as Identity matrices of order 5. Then P=eye(5) shall be the exact value of P which satisfies A2=P*A1*inv(P). How do we obtain this P by solving this equation in matlab code?

[Matt J]

Then P=eye(5) shall be the exact value of P which satisfies A2=P*A1*inv(P).

That will be one solution, but clearly there will be infinitely many alternative solutions as well. Any invertible P would solve the equation.

Answer 1

[Matt J]

A necessary but not sufficient condition for a solution is A2*P-P*A1=0, which is equivalent to K*P(:)=0 where K=kron(I,A2)-kron(A1.',I). The necessary condition can therefore be solved (if a solution exists) with null(). Example:

A1=[1.1581    0.0279   -0.0739   -0.0169   -0.0100
   -0.9172    0.7637    0.5688    0.1299    0.0795
   12.8599    2.5037   -5.4500   -1.4719   -0.8824
  -39.7814   -7.8095   20.0738    5.5809    2.7482
  -27.9052   -5.5369   14.1918    3.2387    2.9446];


P0=rand(5);

A2=P0*A1/P0;

d=length(A1);
I=eye(d);

K=kron(I,A2)-kron(A1.',I);

P=reshape(null(K),d,d,[]);

errors=nan(size(P,3),1);
for j=1:numel(errors)
 Pj=P(:,:,j);
 errors(j)=max(abs(A2-Pj*A1/Pj),[],'all');  
end

P=P(:,:,errors<1e-8)

P = 
P(:,:,1) =

   -0.1484   -0.0809    0.3049    0.1126   -0.0748
   -0.0825    0.4494   -0.5312   -0.1467   -0.2100
   -0.1645    0.1687   -0.0116    0.0257   -0.1725
   -0.1118    0.2680   -0.2137   -0.0442   -0.1774
   -0.1682    0.1182   -0.0494   -0.0055   -0.0969


P(:,:,2) =

    0.1499    0.1345    0.1784    0.0311    0.2248
    0.3233    0.1632   -0.2085    0.2133    0.1304
    0.3195    0.1930   -0.0885    0.0778    0.2600
    0.2801    0.2077   -0.1852    0.1193    0.1836
    0.4238    0.1047   -0.0101    0.0908    0.1241


P(:,:,3) =

    0.2255    0.3731   -0.4906   -0.0203   -0.1859
   -0.2177    0.0766    0.1463    0.1861   -0.1393
    0.0645    0.3138   -0.2742    0.0879   -0.2131
   -0.0570    0.2065   -0.0975    0.1219   -0.1739
    0.0441    0.1757   -0.1544    0.0525   -0.1212


P(:,:,4) =

    0.1178   -0.0384    0.2374    0.1743   -0.2159
    0.2590    0.0815    0.2378    0.1002   -0.2551
    0.2093    0.0079    0.3701    0.2168   -0.3097
    0.2308    0.0235    0.3178    0.1624   -0.2681
    0.0283    0.0140    0.1750    0.0982   -0.1668


P(:,:,5) =

    0.0787   -0.1273   -0.2835   -0.0928   -0.1639
    0.4456   -0.0677   -0.2022   -0.3277   -0.0984
    0.2488   -0.1367   -0.2216   -0.1928   -0.1933
    0.3397   -0.1347   -0.1554   -0.2293   -0.1354
   -0.0409   -0.0693   -0.1722   -0.1530   -0.0908

This does produce several solutions, though it is not clear to me if solving the necessary conditions will always find one.

[Matt J]

This does produce several solutions, though it is not clear to me if solving the necessary conditions will always find one.

Although, if not one could then try fminunc,

Np=null(K);

cols=size(Np,2);

fun=@(coeff)cond(reshape(Np*coeff,d,d));

P=Np*fminunc(fun,rand(cols,1));

[Richa Dubey]

This method isnt working when null(K) in the above code returns an empty matrix. It is working fine for the exaple you showed but not for a general case. Could you please suggest a soln which is 'necessary and sufficient' for satisfying A2=P*A1*inv(P)?

[Matt J]

This method isnt working when null(K) in the above code returns an empty matrix.

If null(K) is empty, then the method has succeeded. It has successfully determined that no solutions exists. For a solution to exist, it must lie in the nullspace of K and therefore null(K) cannot be empty.

Could you please suggest a soln which is 'necessary and sufficient' for satisfying A2=P*A1*inv(P)?

See my comment above about fminunc.

[Richa Dubey]

Thank you so much!

The soln you gave works well but only for those A1 and A2 which satisfy this relation exactly.

Although I was wondering how to deal with those A1 and A2 which don't satisfy A2=P*A1*inv(P) exactly, as in the value of A2 (LHS) and P*A1*inv(P) (RHS) differ by a small error of order say 10^(-3).

How do we calculate P for such A1 and A2 so that this error is minimized. I tried dealing with this using 'lsqnonlin' for a 3rd orded system though but wasn't satisfied with the results.

guess=[1 0.5 0.5 0.5 1 0.5 0.5 0.5 1];
options = optimoptions(@lsqnonlin,'OptimalityTolerance', 1e-100,'Algorithm','levenberg-marquardt', 'FunctionTolerance', 1e-100,'StepTolerance',1e-100,'MaxIterations',1e+100,'MaxFunctionEvaluations',1e+50,'FiniteDifferenceType','central','ScaleProblem','jacobian');
result=lsqnonlin(@eqns,guess,[],[],options);

%result=fsolve(@eqns,guess,options);

P=[result(1)    result(4)   result(7)
    result(2)   result(5)   result(8)
    result(3)   result(6)   result(9)]
function fncs = eqns(z)

    global A1 A2 B1 B2
         a1=z(1); 
         a2=z(2);
         a3=z(3);
         b1=z(4);
         b2=z(5);
         b3=z(6);
         c1=z(7);
         c2=z(8);
         c3=z(9);
         
        P1=[a1 b1 c1
            a2 b2 c2
            a3 b3 c3];
        
       expres=P1*A1*inv(P1)-A2;
       expres1=P1*B1-B2;
       fncs(1)= expres(1,1);
       fncs(2)= expres(1,2);
       fncs(3)= expres(1,3);
       fncs(4)= expres(2,1);
       fncs(5)= expres(2,2);
       fncs(6)= expres(2,3);
       fncs(7)= expres(3,1);
       fncs(8)= expres(3,2);
       fncs(9)= expres(3,3);
       fncs(10)= expres1(1);
       fncs(11)= expres1(2);
       fncs(12)= expres1(3);
       
end

I calculated error using this:

A22=P*A1*inv(P);
B22=P*B1;
eA=((A22-A2)./A2)*100
eB=((B22-B2)./B2)*100

Answer 2

[Torsten]

I think you could also get such P by solving the linear optimization problem

min: sum_ij (eij+ + eij-)

under the constraints

E+ - E- = (P+ - P-)*A1 - A2*(P+ - P-)
sum_ij (pij+ + pij-) = 1
E+,E-,P+,P- >= 0

The matrix P = P+ - P- should have the property P*A1 = A2*P if the optimal value of the objective is 0.

Otherwise, such P does not exist.

Note that it's not guaranteed that P is invertible (I guess). But it's different from the trivial solution P=0.

[Richa Dubey]

Thanks for the soln.

I am unable to implement it though. Could you please elaborate more on this?

Also, please look into my comment in the above ans if possible.