# 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. – Henrik Schumacher Aug 26 at 11:51
• The solution to $A = B^{-1}$ seems to be $B=A^{-1}$, no? – Michael E2 Aug 26 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$. – Mathieu Aug 26 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}$. – Mathieu Aug 26 at 12:15

## 1 Answer

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.

• Thanks :) As I am new to Mathematica, could you detail a bit what the code does? Are you also deleting the values below the diagonal? (which is actually great because those equations are redundant). – Mathieu Aug 26 at 12:49
• What is the asset of transforming the equations into polynomial instead of the rational? – Mathieu Aug 26 at 12:55
• 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? – Mathieu Aug 26 at 14:32
• en.wikipedia.org/wiki/Cramer%27s_rule – Henrik Schumacher Aug 26 at 14:34
• 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/… – Mathieu Aug 27 at 9:13