# 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?

1 view (last 30 days)
Richa Dubey on 22 Jan 2022
Commented: Richa Dubey on 23 Jan 2022
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!
##### 2 CommentsShowHide 1 older comment
Matt J on 22 Jan 2022
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.

Matt J on 22 Jan 2022
Edited: Matt J on 22 Jan 2022
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.
Richa Dubey on 23 Jan 2022
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

Torsten on 22 Jan 2022
Edited: Torsten on 22 Jan 2022
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 on 23 Jan 2022
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.

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!