# How to speed up filling a matrix as much as possible in a loop?

In order to find the solution of an $n\times n$ system of nonlinear equations by the method $x_{k+1}=x_k-D_k^{-1}f(x_k)$ where $D_k$ is an approximation of the Jacobian and defined by

\begin{equation} [x,y;F]_{i,j}=\frac{F_i(x_1,\ldots,x_j,y_{j+1},\ldots,y_n)-F_i(x_1,\ldots,x_{j-1},y_{j},\ldots,y_n)}{x_j-y_j}, \quad 1\leq i,j\leq n, \end{equation} this matrix must be defined in each loop. Although I defined it below for an example in high precision, it seems very slow!

My code is in what follows:

 (*Defining the nonlinear system and initial guess*)
ClearAll["Global*"];
digits = 400; size = 10;
Id = SparseArray[{{i_, i_} -> 1}, {size, size}];
x = SetAccuracy[{1.4 + 0.5 I, 1. - 2.0 I, 1.0 - 0.2 I, 2.5 + 0.5 I,
0.8 - 0.1 I, -0.4 + 1. I, 0.1 + 0.1 I, 1.3 - 0.7 I, 2.0 + 0.5 I, -2.0 + 1.5 I},
digits];
f[{x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_, x10_}] :=
{5 Exp[x1 - 2]*x2 + 8 x3^x4 - 5 x6^3 + 2*x7^x10 - x9,
5 Tan[x1 + 2] + x2^3 + 7 x3^4 - 2 Sin[x6]^3 + Cos[x9^x10],
x1^2 + Tan[x2] + 2 x3^x4 - 5 x6^3 - x5*x6*x7*x8*x9*x10,
2 Tan[x1^2] + 2^x2 + x3^2 - 5 x5^3 - x6 + x8^Cos[x9],
10 x1^2 + Cos[x2] + x3^2 - 5 x6^3 - 4^x9 - 2 x8 - x10,
ArcCos[x1^2]*Sin[x2] + x3^2 - 2 x5^4*x6*x9*x10,
x1*x2^x7 + x3^5 - 5 x5^3 + x7 - x8^x10,
x4*Sin[x2] + x3 - 15 x5^2 + x7 + ArcCos[x8 + x9 - 10 x10],
10 x1 + x3^2 - 5 x5^2 + 10 x6^x8 + 2 x9 - Sin[x7],
x1*Sin[x2] - 5 x6 - 2 x10^x8 - 10 x9 + x10};


And also

 (*Defining the iteration loop*)
point1[j_] := Flatten@{Table[x[[k]], {k, 1, j}], Table[w[[k]], {k, j + 1, size}]};
point2[j_] := Flatten@{Table[x[[k]], {k, 1, j - 1}], Table[w[[k]], {k, j, size}]};
max = 11; fx = f[x];
β =  SparseArray[{{i_, i_} -> -1/100}, {size, size}];
Do[{W = w; w = x + β.fx; fw = f[w];
T = Table[(f[point1[j]][[i]] - f[point2[j]][[i]])/(x[[j]] - w[[j]]),
{i, size}, {j, size}];
yy = LinearSolve[T, fx];
X = x; x = x - yy; fX = fx; fx = f[x];
L[i] = Norm[fx, 2];},
{i, 1, max}]; // AbsoluteTiming

Table[N[ScientificForm[L[k]], 4], {k, 1, max}]
N[x, 5]


Unfortunately, this process takes around 7 seconds in my laptop and I believe if I could fill in the matrix $T$ in the code faster, then the whole time will be reduced.

Note that $\mathtt{ParallelTable[]}$ reduces the computational time but not significantly!

I will be thankful if anyone could give me some tips of how to fill in the matrix as quickly as possible.

The following is ten times faster. The remaining time is mostly consumed while evaluating your function, so there may be some optimization window there.

point1[j_] := Join[x[[;; j]], w[[j + 1 ;;]]];
point2[j_] := Join[x[[;; j - 1]], w[[j ;;]]];
max = 11;
fx = f[x];
β = SparseArray[{{i_, i_} -> -1/100}, {size, size}];
Do[{
w = x + β.fx;
T = Transpose@Table[(f[point1[j]] - f[point2[j]])/(x[[j]] - w[[j]]), {j, size}];
x = x - LinearSolve[T, fx];
fx = f[x];
L[i] = Norm[fx, 2];}, {i, 1, max}]; // AbsoluteTiming

Table[N[ScientificForm[L[k]], 4], {k, 1, max}]
N[x, 5]
`
• Wow, your answer is unbelievable. It works perfectly and speed up the process of filling the matrix. Thanks a lot. – Fazlollah Jan 14 '16 at 19:21
• I have a related observation. For large systems ($n>200$), once again filling the matrix is slow. Do you have any ideas to fill the matrix for sparse cases also too fast? – Fazlollah Jan 19 '16 at 19:55
• @FazlollahSoleymani Are you sure that the cost of calculating the function isn't the dominant one? – Dr. belisarius Jan 19 '16 at 19:58
• Yes, the cost of calculating the function is the dominant one. But, when the size of the system becomes large, computing the function evaluations and then filling the matrix becomes slow again. I am trying to solve a discretization nonlinear system resulting from a PDE with the above approach. But for large sizes, it is slow. Maybe, it would be better to incorporate the notion of sparsity in your nice answer. – Fazlollah Jan 20 '16 at 14:25
• @fazl in what sense is the system "sparse"? – Dr. belisarius Jan 21 '16 at 2:14