-2

how can I solve this overdetermined equation system with the command PseudoInverse. enter image description here

left = Table[(-i)^i, {i, 1, 4}]; right1 = Table[Sum[2^(i*k)*Sin[Subscript[x, k]], {k, 1, 3}], {i, 1, 4}];

 right2 =Table[Sum[2^(i*k)*Cos[Subscript[x, k]], {k, 2, 3}], {i, 1, 4}]; 
 StringForm["left side = ``", left // MatrixForm] 
 StringForm["right1 side = ``", right1 // MatrixForm]
 StringForm["right2 side = ``", right2 // MatrixForm]
 A1 =  {{2, 4, 8},{4, 16, 64},{8, 64, 512},{16, 256, 4096}};  
 A2 = {{4, 8},{16, 64},{64, 512},{256, 4096}};
 b ={{-1},{4},{-27},{256}};
 x1 = PseudoInverse[Subscript[A, 1]].b // N;
 x2 = PseudoInverse[Subscript[A, 2]].b // N;
 r1 = Subscript[A, 1].Subscript[x, 1] - b // MatrixForm;
 r2 = Subscript[A, 2].Subscript[x, 2] - b // MatrixForm;
 r1 = Flatten[{{r1[[1, 1]]},{r1[[1, 2]]},{r1[[1, 3]]},{r1[[1, 4]]}}];
 r2 = Flatten[{r2[[1, 1]]},{r2[[1, 2]]},{r2[[1, 3]]},{r2[[1, 4]]}}];
 g1 = Length[Subscript[r, 1]];
 g2 = Length[Subscript[r, 2]];
 p1 = Sqrt[Subscript[r, 1].Subscript[r, 1]/Subscript[g, 1]];
 p2 = Sqrt[Subscript[r, 2].Subscript[r, 2]/Subscript[g, 2]];
 Print["r1 = ", r1 // MatrixForm, "\np1 = ", p1]
 Print["r2 = ", r2 // MatrixForm, "\np2 = ", p2]
  • 1
    Please, share the code in copyable and syntactically correct form. Also, this is not a linear system of equations, so Pseudoinverse cannot be applied directly. The method of choice would proabably be the Gauß-Newton algorithm (or rather the Levenberg-Marquardt variant). It is employed by FindMininum if you reformulate your problem as a regression or parameter fitting problem. FindFit and friends might also be of help here. – Henrik Schumacher Apr 05 '19 at 08:34

1 Answers1

4

This is not a linear system of equations, so Pseudoinverse cannot be applied directly. The method of choice would proabably be the Gauß-Newton algorithm (or rather its Levenberg-Marquardt variant). Under the hood, it employs pseudoinverses of the linearizations of the system.

The Levenberg-Marquardt method is employed by FindMininum and probably also by NMinimize if one reformulates the problem as a regression or parameter fitting problem:

eq = Table[(-i)^(-i) - Sum[(-2)^(k i) If[k == 1, Cos[x[k]], Sin[x[k]]], {k, 1, 3}], {i, 1, 4}];
sol = NMinimize[Total[eq^2], Table[x[k], {k, 1, 3}]]

{0.709009, {x[1] -> 1.44635, x[2] -> 0.00659458, x[3] -> -0.000899084}}

As you can see, a solution to your system could not be found, as the minimal value is 0.709009 rather far away from 0..

FindFit and friends might also be of help here.

Edit

Daniel's comment below rang a bell. The following show that not solution exists because the right hand side of the equations does not lie within the image of the matrix A (because it is not perpendicular to the orthogonal complement of the image of A which is the null space of Transpose[A]):

vars = {Cos[x[1]], Sin[x[2]], Sin[x[3]]};
left = Table[Sum[(-2)^(k i) vars[[k]], {k, 1, 3}], {i, 1, 4}];
right = Table[(-i)^(-i), {i, 1, 4}];
A = D[left, {vars, 1}];

NullSpace[A[Transpose]][[1]].right

133129/2304

Actually,

Solve[left == right, {x[1], x[2], x[3]}]

{}

tells us the same.

If the right hand side were in the image of A, all solutions of the equations could be obtained by the following:

pinv = PseudoInverse[A].right;
Solve[vars == pinv, {x[1], x[2], x[3]}]

Notice that Nullspace[A] returns {}, so A is injective and these would be indeed all solutions. In our case, there are no solutions of your equations, so what Solve returns here is the set of least squares solutions.

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309