I am desperately trying to resolve an overdetermined system using Maxima and the lsquares_estimate function.
The system is described by this matrix equation : M.D = I
Here is M :
M:matrix([19.707,0.06700000000000017,0,0,0,0,0,0,19.707-Phi01_8,0,0,0],
[0,0,-0.2690000000000019,0.7359999999999971,0,0,0,0,0,19.371-Phi01_8,0,0],
[-0.067234,-24.562234,0,0,0,0,0,0,-Phi12_8-0.067234,0,0,0],
[0,0,-1.753,3.741,0,0,0,0,0,22.742-Phi12_8,0,0],
[-0.26894,1.75226,0,0,0,0,0,0,-Phi23_8-0.26894,0,0,0],
[0,0,92.99119999999999,14.967,0,0,0,0,0,90.97-Phi23_8,0,0],
[-0.73571,3.74149,0,0,0,0,0,0,-Phi34_8-0.73571,0,0,0],
[0,0,-14.9658,-655.2429999999999,0,0,0,0,0,-Phi34_8-19.443,0,0],
[-1.8393,9.3537,0,0,0,0,0,0,-Phi45_8-1.8393,0,0,0],
[0,0,-37.415,361.892,0,0,0,0,0,-Phi45_8-48.608,0,0],
[-8.6658,5.6092,0,0,0,0,0,0,-Phi56_8-8.6658,0,0,0],
[0,0,-22.437,166.598,0,0,0,0,0,-Phi56_8-36.712,0,0],
[-5.1995,3.3655,0,0,0,0,0,0,-Phi67_8-5.1995,0,0,0],
[0,0,-13.462,99.953,0,0,0,0,0,-Phi67_8-22.027,0,0])
Here is D :
D:matrix([0.0499998750003125],
[0.0399993670651863],
[0.009999946119925522],
[0.001250003764581481],
[4.999947987060861*10^-4],
[1.000099826698366*10^-4],
[1.666560996027836*10^-4],
[0.00499998750003125],
[D8],[D9],[D10],[D11])
And here is I :
I:matrix([1],[0],[-1],[0],[0],[1],[0],[-1],[0],[0],[0],[0],[0],[0])
The system has 14 equations and 9 variables. (D8,D9,Phi01,Phi12,Phi23,Phi34,Phi45,Phi56,Phi67) I want to get a least squares estimate of those variables. I am using the "lsquares_estimates" function. In order to give the right input to the function, I merge the matrix equation into one matrix. This is how i call the lsquares_estimates function :
D0:D[1][1]$
D1:D[2][1]$
D2:D[3][1]$
D3:D[4][1]$
D4:D[5][1]$
D5:D[6][1]$
D6:D[7][1]$
D7:D[8][1]$
D8:D[9][1]$
D9:D[10][1]$
D10:D[11][1]$
D11:D[12][1]$
C:addcol(M,I);
resultat:lsquares_estimates(C,
[c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,Ie],
Ie=c0*D0+c1*D1+c2*D2+c3*D3+c4*D4+c5*D5+c6*D6+c7*D7+c8*D8+c9*D9+c10*D10+c11*D11,
[Phi01_8,Phi12_8,Phi23_8,Phi34_8,Phi45_8,Phi56_8,Phi67_8,D8,D9]);
Obviously, Maxima is not able to resolve this system... It keeps running for a very longtime, and sometimes, here is the error i get :
Maxima encountered a Lisp error:
-: ((MTIMES SIMP) -0.14285714285714285
((MPLUS SIMP)
((MTIMES SIMP)
((MPLUS SIMP) 0.01197250577547404
((MTIMES SIMP) -1.0 ((MPLUS SIMP) 19.707 ((MTIMES SIMP) -1.0 |$Phi01_8|))))
((MPLUS SIMP) 19.707 ((MTIMES SIMP) -1.0 |$Phi01_8|)))
((MTIMES SIMP)
((MPLUS SIMP) -0.014164494697229801
((MTIMES SIMP) -1.0
((MPLUS SIMP) -0.067234 ((MTIMES SIMP) -1.0 |$Phi12_8|))))
((MPLUS SIMP) -0.067234 ((MTIMES SIMP) -1.0 |$Phi12_8|)))
((MTIMES SIMP)
((MPLUS SIMP) -0.056642324551059296
((MTIMES SIMP) -1.0
((MPLUS SIMP) -0.26894 ((MTIMES SIMP) -1.0 |$Phi23_8|))))
((MPLUS SIMP) -0.26894 ((MTIMES SIMP) -1.0 |$Phi23_8|)))
((MTIMES SIMP)
((MPLUS SIMP) -0.11287182384424398
((MTIMES SIMP) -1.0
((MPLUS SIMP) -0.73571 ((MTIMES SIMP) -1.0 |$Phi34_8|))))
((MPLUS SIMP) -0.73571 ((MTIMES SIMP) -1.0 |$Phi34_8|)))
((MTIMES SIMP)
((MPLUS SIMP) -0.28217730962955834
((MTIMES SIMP) -1.0
((MPLUS SIMP) -1.8393 ((MTIMES SIMP) -1.0 |$Phi45_8|))))
((MPLUS SIMP) -1.8393 ((MTIMES SIMP) -1.0 |$Phi45_8|)))
((MTIMES SIMP)
((MPLUS SIMP) 0.2089244670356651
((MTIMES SIMP) -1.0
((MPLUS SIMP) -8.6658 ((MTIMES SIMP) -1.0 |$Phi56_8|))))
((MPLUS SIMP) -8.6658 ((MTIMES SIMP) -1.0 |$Phi56_8|)))
((MTIMES SIMP)
((MPLUS SIMP) 0.12535648020624035
((MTIMES SIMP) -1.0
((MPLUS SIMP) -5.1995 ((MTIMES SIMP) -1.0 |$Phi67_8|))))
((MPLUS SIMP) -5.1995 ((MTIMES SIMP) -1.0 |$Phi67_8|))))) is not a number
Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
Am I using lsquares_estimates the wrong way ? When i reduce the number of equations to match the number of variables, I can get a result, but i would like an accurate result using all my equations.
Any help would be really appreciated... Thanks in advance.
lsquares_estimates constructs a mean square error expression and then tries to solve it exactly, and then approximately, if solving it exactly doesn't yield a result. Not sure where things are going wrong there.
But anyway you can jump directly to lsquares_estimates_approximate. Try this:
e1 : M.D - I;
my_mse : transpose(e1) . e1;
load (lsquares);
lsquares_estimates_approximate (my_mse, [Phi01_8, Phi12_8, Phi23_8,
Phi34_8, Phi45_8, Phi56_8, Phi67_8, D8, D9], tol=1e-4);
Note that my_mse is the squared norm of the error.
With these inputs I get:
[[Phi01_8 = 22.70059866292853,Phi12_8 = -4.329661712144546,
Phi23_8 = -17.40507905561839,Phi34_8 = -34.32504757490423,
Phi45_8 = -85.85849926978096,Phi56_8 = 54.09091751735507,
Phi67_8 = 32.46719922407783,D8 = -0.003346473781793617,
D9 = 2.1359057759474724E-5]]