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

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?

• Impossible to say without the any detailed code. Aug 26 '19 at 11:51
• The solution to $A = B^{-1}$ seems to be $B=A^{-1}$, no? Aug 26 '19 at 12:03
• @MichaelE2 Ok, it seems I poorly explained the problem in this post. Sorry. $A$ is a matrix where there are no $0$ at the beginning. But, before the solving, I replace the values $A_{i,j}$ with $0$ IF $B_{i,j} = 0$. Then, $B^{-1}$ elements are all a non linear combination of the $b_{i,j}$, the unkowns that I am looking for. So the solving as too solve element-wise the problem to figure out what the $b_{i,j}$ are. In other words, if I put $0$ in $A$ where $B_{i,j}=0$, and then invert $A$, I won't get a $0$ anymore in the location $A^{-1}_{i,j}$ where $B_{i,j}=0$. Aug 26 '19 at 12:11
• @HenrikSchumacher I added the matrix definition, but all the code is here. I am simply defining $A$ and $B$, then I invert $B$ which gives a matrix of equations, and finally, I create 2 tables from $A$ and $Binv$ where a $0$ is placed at the location $(i, j)$ if $B_{i,j} = 0$. Then I solve the equality between those tables. As NSolve is not converging, I would like to use FindRoot where I define the starting point for the variable $b_{i,j}$ as $A^{-1}_{i,j}$. Aug 26 '19 at 12:15

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];
• One more question :) Why is the fraction obtained from the inversion $eq / det(B)$? Is this due to the form of the matrix $B$? Will it always be the case no matter where I have $B_{i,j} = 0$ as long as $B$ is symmetric and $B_{i,i} /neq 0$. i.e. will it always be the $det(B)$ if I have symmetric $B$ with a filled diagonal and some $0$ on the off-diagonal elements? Aug 26 '19 at 14:32
• Thanks for all the explanation. While going full scale I ran into a problem. The computation of the determinant of $B$ (16x16) is taking forever. Is this step really needed since you say that it doesn't make a difference in this case? I added the 16x16 matrix in this post: mathematica.stackexchange.com/questions/204374/… Aug 27 '19 at 9:13