Numerically solve system of 64 equations with 64 variables in matrix form

Question

I have 2 matrix $A$ and $B$ both of size 16x16. Most of the values $A_{i,j}$ and $B_{i,j}$ for the same couple $(i,j)$ are equal to 0. The remaining ones are a positive floating number for $A$ and a non-linear combination of the 64 variables for $B$.

In other words, for any $(i,j)$, $A_{i,j} = B_{i,j} = 0$ OR $A_{i,j} = number$ and $B_{i,j} =$ $non-linear$ $combination$ $of$ $the$ $64$ $variables$.

In the end, I have a 64 equations with 64 variables equations system in the form $A = B$.

I am trying to solve this with NSolve on a workstation with 128 Gb of RAM. After 72 hours of computation at nearly 120 Gb of RAM non-stop... I'm starting to wonder if it will ever converge...

I just read that another function may be performing better FindRoot. I do have an idea of the value for each of the 64 variables, however, I do not know how I can input those starting points for the 64 variables.

4x4 example:

$$A = \begin{pmatrix}50 & 3 & 10 & 2\\\ 3 & 60 & 7 & 1\\\ 10 & 7 & 55 & 4\\\ 2 & 1 & 4 & 45 \end{pmatrix}$$

$$B = \begin{pmatrix}b_{11} & b_{12} & b_{13} & 0\\\ b_{21} & b_{22} & 0 & 0\\\ b_{31} & 0 & b_{33} & b_{34}\\\ 0 & 0 & b_{43} & b_{44} \end{pmatrix}$$

A = {{50, 3, 10, 2}, {3, 60, 7, 1}, {10, 7, 55, 4}, {2, 1, 4, 45}}
B = {{Subscript[b, 11], Subscript[b, 12], Subscript[b, 13], 
  0}, {Subscript[b, 12], Subscript[b, 22], 0, 0}, {Subscript[b, 13], 
  0, Subscript[b, 33], Subscript[b, 34]}, {0, 0, Subscript[b, 34], 
  Subscript[b, 44]}}

Solve $A = B^{-1}$.

Binv = Inverse[B]
NSolve[Table[
   If[MatchQ[B[[i, j]], Subscript[b, x_]], Binv[[i, j]], 0], {i, 
    4}, {j, 4}] == 
  Table[If[MatchQ[B[[i, j]], Subscript[b, x_]], A[[i, j]], 0], {i, 
    4}, {j, 4}]]

Solution:

$$B=\left(\begin{array}{cccc}\frac{27497}{1321025} & -\frac{1}{997} & -\frac{1}{265} & 0 \\-\frac{1}{997} & \frac{50}{2991} & 0 & 0 \\-\frac{1}{265} & 0 & \frac{136093}{7167985} & -\frac{4}{2459} \\0 & 0 & -\frac{4}{2459} & \frac{55}{2459} \\\end{array}\right)$$

I would like to find the same solution with FindRoot by using the values $A^{-1}_{i,j}$ as starting point for $b_{i,j}$ but I don't know how to provide this input to FindRoot.

Additionnaly, are there also any other parameters I should provide to FindRoot?

Answer 1

At least in your toy model, using the respective entries of the inverse of A for initializing FindRoot seems to work quite fine:

B0 = B /. Subscript[b, _] :> 1;
eq = DeleteCases[Flatten[B0 Det[B] UpperTriangularize[Inverse[B] - A]], 0];
Binit = DeleteCases[Flatten[UpperTriangularize[B0 Inverse[A]]], 0];
vars = DeleteCases[Flatten[UpperTriangularize[B]], 0];
FindRoot[Thread[eq == 0], Transpose[{vars, Binit}]]

{Subscript[b, 11] -> 0.0208149, Subscript[b, 12] -> -0.00100301, Subscript[b, 13] -> -0.00377358, Subscript[b, 22] -> 0.0167168, Subscript[b, 33] -> 0.0189862, Subscript[b, 34] -> -0.00162668, Subscript[b, 44] -> 0.0223668}

I remove the redundant equations with UpperTriangularize and DeleteCases because otherwise FindRoot would complain about the fact that there are more equations than variables.

I also multiply the equations with Det[B]; this way the equations become polynomial -- beforehand they were rational. Rational equations can be problematic for Newton's method (which is employed by FindRoot) because it can lead to division by zero. In this case, it doesn't make a difference, though.