# How to speed up auxilary DoolittleDecomposite function?

The Doolittle Decomposition algorithm as below: $$A=LU$$ where

$A= \begin{pmatrix} a_{11} & a_{12} & \cdots a_{1n}\\ a_{21} & a_{22} & \cdots a_{2n}\\ \vdots & \vdots & \vdots \\ a_{n1} & a_{n2} & \cdots a_{nn}\\ \end{pmatrix} L=\begin{pmatrix} 1 & 0 & \cdots 0\\ l_{21} & 1 & \cdots 0\\ \vdots & \vdots & \vdots \\ l_{n1} & l_{n2} & \cdots 1\\ \end{pmatrix}$ $U=\begin{pmatrix} u_{11} & u_{12} & \cdots u_{1n}\\ 0 & u_{12} & \cdots u_{12}\\ \vdots & \vdots & \vdots \\ 0 & 0 & \cdots u_{nn}\\ \end{pmatrix}$

And

\begin{equation} \begin{cases} u_{kj}=a_{kj}-\sum_{m=1}^{k-1}l_{km}u_{mj} && (j=k,k+1,\cdots,n)\\ l_{ik}=1/u_{kk} (a_{ik}-\sum_{m=1}^{k-1}l_{im}u_{mk}) && (i=k+1,k+2,\cdots,n) \end{cases} \end{equation}

Where, $k=1,2,\cdots,n$

### Implementation

oneStep[mat_?MatrixQ, {intermediate_, index_}] :=
Module[{temp = intermediate, k = index, row},
row = Length@mat;
Table[
temp[[k, j]] =
mat[[k, j]] - Sum[temp[[k, m]] temp[[m, j]], {m, 1, k - 1}], {j, k, row}];
Table[
temp[[i, k]] =
1/temp[[k, k]] (mat[[i, k]] -
Sum[temp[[i, m]] temp[[m, k]], {m, 1, k - 1}]), {i, k + 1, row}];
{temp, k+1}
]
(*==============================*)

doolittleDecomposite[mat_?MatrixQ] :=
Module[{intermediate, row},
intermediate = ConstantArray[0, Dimensions@mat];
row = Length@mat;
First@
Nest[oneStep[mat, #] &, {intermediate, 1}, row]
]


### Test

test={{1, 6, 1, 9}, {4, 4, 9, 7}, {10, 4, 10, 2}, {10, 3, 10, 2}};
doolittleDecomposite[test]


$\left( \begin{array}{cccc} 1 & 6 & 1 & 9 \\ 4 & -20 & 5 & -29 \\ 10 & \frac{14}{5} & -14 & -\frac{34}{5} \\ 10 & \frac{57}{20} & \frac{57}{56} & \frac{11}{7} \\ \end{array} \right)$

Namely,

$L=\left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 4 & 1 & 0 & 0 \\ 10 & \frac{14}{5} & 1 & 0 \\ 10 & \frac{57}{20} & \frac{57}{56} & 1 \\ \end{array} \right) U=\left( \begin{array}{cccc} 1 & 6 & 1 & 9 \\ 0 & -20 & 5 & -29 \\ 0 & 0 & -14 & -\frac{34}{5} \\ 0 & 0 & 0 & \frac{11}{7} \\ \end{array} \right)$

### Performance testing

SeedRandom;
testReal = RandomReal[{1., 100.}, {200, 200}];
doolittleDecomposite[testReal]; // AbsoluteTiming
LUDecomposition[testReal]; // AbsoluteTiming

{6.102349, Null}
{0.034002, Null}  (*big performance difference*)

SeedRandom
testInteger = RandomInteger[{1, 100}, {100, 100}];
doolittleDecomposite[testInteger]; // AbsoluteTiming
LUDecomposition[testInteger]; // AbsoluteTiming

{6.173353, Null}
{4.481256, Null}  (*a bit difference*)


### Question

• Is it possible to speed up my implementation code?
• Could someone give me some hints(like how to organize the data) to implement this algorithm elegantly in Mathematica?Namely, I think my method is not good:)
• Might try changing the sums to dot products. Also Table probably should be Do. For the machine real case you could run it through Compile. – Daniel Lichtblau Mar 31 '15 at 16:05
• The oneStep function that I implemented is C-style(like procedural paradigm) . So my main difficult is how to store and organize the intermediate result in right way if I'd like to implement this algorithm with functional or rule-based paradigm. – xyz Apr 3 '15 at 5:21
• Might want to check the notebook here. Dated a bit but still quite relevant to the question at hand. – Daniel Lichtblau Apr 4 '15 at 13:20
• I've just seen this. I believe you know, but just to remind, algorithms of this sort are usually made with a pivoting strategy in place (e.g. partial, total, or rook); otherwise, there are a lot of nonsingular matrices that do not have a decomposition. – J. M. is in limbo Oct 24 '15 at 10:56
• With respect to the counterexample, try Reverse[IdentityMatrix]. As for pivoting: did you notice the difference between your method and LUDecomposition[]? The latter one outputs a permutation list, which encodes the (partial) pivoting (interchanges of rows) done in the course of decomposing the matrix. Any robust implementation has to have a pivoting strategy like that. – J. M. is in limbo Oct 24 '15 at 11:55

It's pretty clear that the complexity of this function is $\text{O}(n^3)$, since we're iterating over each element of a matrix (a factor of $n^2$) and taking a sum at each one (an additional factor of $n$). As a former USACO competitor, I don't like $\text{O}(n^3)$ algorithms. However, we can't really do anything to reduce the complexity; LU decomposition has the same complexity as matrix multiplication. We can get faster than $\text{O}(n^3)$: based on my testing, Mathematica uses an $\text{O}(n^{\log_2 7})$ multiplication algorithm, probably this one; but such algorithms are complex, memory-intensive, and beyond the scope of this question.

The upshot of all this is that if we really want to increase our speed, we should Compile our function. Here is an algorithm in Fortran that is in-place; that is, it doesn't use any extra memory (although we do need to copy the matrix once so that we can 'edit' it). We can write a CompiledFunction to implement this version of the algorithm:

doolittleDecomposition = Compile[{{a, _Real, 2}},
With[{n = Length[a]},
Module[{t = a},
Do[
Do[
Do[
t[[i, j]] -= t[[i, k]] t[[k, j]];,
{k, 1, i - 1}];,
{j, i, n}];
Do[
Do[
t[[j, i]] -= t[[j, k]] t[[k, i]];,
{k, 1, i - 1}];
t[[j, i]] /= t[[i, i]];,
{j, i + 1, n}];,
{i, 1, n}];
t
]
],
{{q, _Real}}
, RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False},
CompilationTarget -> "C",
CompilationOptions -> {"ExpressionOptimization" -> True}]


We should double-check that our function doesn't call MainEvaluate, which is slow:

Needs["CompiledFunctionTools"]
doolittleDecomposition // CompilePrint


While doing this check I discovered that I originally wrote t[j, i] instead of t[[j, i]], causing MainEvaluate to be called. Always double-check! In fact, let's now double-check that both algorithms give the same result:

testReal = RandomReal[1, {200, 200}];
d1 = doolittleDecomposite[testReal]; // AbsoluteTiming
d2 = doolittleDecomposition[testReal]; // AbsoluteTiming
d1 == d2


This also shows that the Compiled version is faster.

## Update: Refactoring the original code

If you think Compile is too drastic (or you don't have a working C compiler) there are some other changes we can make.

• Don't use an auxiliary or 'helper' function. This will incur extra overhead that you don't need. Unlike a compiled language where the compiler will essentially inline the function, Mathematica will be passing-by-value all of your matrices!
• Don't use Nest. If we're getting rid of the helper function, Nest doesn't really make sense. We can replace it with a Do-loop instead. The Do-loop can 'automatically' increment k, instead of having to pass it (in a list!) to a function and have it returned (also in a list!!). At this point we can rename intermediate to temp, consolidating two variables that really contain the same thing, and we also can lose the extra redefinition of row inside the loop.
• Use Dot instead of Sum. In general, if there's a builtin function for something, that function is faster. In this case you want to pairwise multiply two lists and take the sum of the resulting elements. We usually call that a dot product. Instead of looping m, we can use Span to extract a list using Part, and then Dot our two lists together. This is where most of the speedup comes from.

At this point the code looks pretty much like the Fortran function, but with the inner loops handled by Dot:

doolittleDecomposite2[mat_?MatrixQ] :=
Module[{temp = ConstantArray[0, Dimensions@mat], row = Length@mat},
Do[
Do[temp[[k, j]] =
mat[[k, j]] - temp[[k, ;; k - 1]].temp[[;; k - 1, j]],
{j, k, row}];
Do[temp[[i, k]] =
(mat[[i, k]] - temp[[i, ;; k - 1]].temp[[;; k - 1, k]]) / temp[[k, k]],
{i, k + 1, row}];,
{k, row}];
temp
]


## Updated performance analysis

I ran another comparison, this time including LUDecomposition, just to see how well we did: The first thing I noticed (even before I generated the plots) was that the original solution takes an unreasonably long time on matrices larger than around $n>250$ (I tried $n=400$ and let it run for at least an hour before aborting manually). I suspect this has to do with memory bandwidth when passing temp into Sum.

The second think to notice is that LUDecomposition is fast. Not only is it at least two orders of magnitude faster than mine, but below $n=1000$ it is almost instant (less than 10 ms). The former is probably from using the Intel Math Kernel Library, which, unlike our code, will max out all 8 threads on my CPU, instead of running single-threaded. It's also optimized to fully utilize every CPU cycle by utilizing advanced instructions and avoiding cache misses and pipeline stalls. As for the latter, it's probably not a coincidence that the 'fast' runs all have memory usage below the size of my CPU cache: The third thing to notice is that while all of our solutions are firmly $\text{O}(n^3)$ for large numbers, the native LUDecomposition is actually closer to $\text{O}(n^{2.4})$, meaning that it must be based on some variant of the Coppersmith–Winograd algorithm (the fastest known matrix-multiplication method). (one of the developers said this is not true) there is probably a large $\text{O}(n^2)$ factor that the leading $n^3$ won't dominate until very large $n$.

Although the refactored version appears to be slightly better than $\text{O}(n^3)$, this is only because of all the overhead; as the vectors passed to Dot get larger, the function becomes asymptotically more efficient, and the $n^3$ term dominates.

Looking back at the memory usage plot, we can see that the original function doesn't use much more memory than mine. All are asymptotically $\text{O}(n^2)$. The reason the original takes so long is because Mathematica is essentially pass-by-value, so it's making a huge malloc and free every time it calls the helper function.

The refactored version uses almost exactly the same amount of memory as my Compiled solution (excepting a small $\text{O}(n)$ contribution from the vectors and Dot) because we avoid all that argument-passing and copying.

As to why I'm using almost exactly 3× more memory than LUDecomposition, it's probably also due to pass-by-value semantics: Mathematica constructs a version of the matrix to pass to the CompiledFunction, which then makes it's own working copy, which is then copied again to produce the result. As LUDecomposition is a native function, it can avoid all that.

• @ShutaoTang does the update answer your question about "how to store and organize the intermediate result?" (Namely, don't copy it around.) – 2012rcampion Apr 4 '15 at 4:15
• You can use CompileGetElement to speed up doolittleDecomposition by a factor of 2. (See here for an example. ) – xzczd Apr 4 '15 at 4:25
• +1Nice answer,thanks a lot. In fact, I'd like to implement this algorithm in functional paradigm(like Nest, because I think functional paradigm is more efficient) at the beginning. However, I didn't know how to organize the intermediate result(like $u_{11},u_{12} \cdots u_{1n}$). So I use the Do(I am always thinking it is not efficient in Mathematica) to implement it. In addition, @Daniel Lichtblau gave the suggestion(changing the Sum to dot products), but I forgot to use the built-in functionSpan(;;) . So your edit fo my original code is very elegant:) – xyz Apr 4 '15 at 4:35
• @xzczd Interesting, why is that not the default behavior of Compile? – 2012rcampion Apr 4 '15 at 4:35
• @ShutaoTang The Doolittle algorithm is essentially procedural: shoehorning it into a functional construct will not make it faster. The reason we advise a functional style most of the time is that builtins are usually faster; but there are no builtins that apply to this case. I believe my testing shows that here, Do is faster than Nest. – 2012rcampion Apr 4 '15 at 4:40

The function doolittleDecomposite2 refactored by @2012rcampion can use Span(;;) to avoid the inner Do loop

doolittle[mat_?MatrixQ] :=
Module[
{temp = ConstantArray[1, Dimensions@mat], row = Length@mat},
Do[
temp[[k, k ;; row]] =
mat[[k, k ;; row]] - temp[[k, ;; k - 1]].temp[[;; k - 1, k ;; row]];
temp[[k + 1 ;; row, k]] =
(mat[[k + 1 ;; row, k]] -
temp[[k + 1 ;; row, ;; k - 1]].temp[[;; k - 1, k]])/temp[[k, k]];,
{k, row}];
temp
]


Or

LU[mat_] :=
Module[{n, L, U},
n = Length@mat;
L = IdentityMatrix[n];
U = mat;
Do[
L[[k ;; n, k]] = U[[k ;; n, k]]/U[[k, k]];
U[[(k + 1) ;; n, k ;; n]] =
U[[(k + 1) ;; n, k ;; n]] - L[[(k + 1) ;; n, {k}]].U[[{k}, k ;; n]];
, {k, 1, n - 1}];
{L, U}
]

(*=========================*)
mat = RandomInteger[10, {5, 5}]
LUDecomposition@mat


$\left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \\ \frac{4}{7} & 1 & 0 & 0 & 0 \\ 1 & -\frac{14}{5} & 1 & 0 & 0 \\ 1 & -\frac{49}{20} & \frac{219}{176} & 1 & 0 \\ \frac{5}{7} & -\frac{1}{2} & \frac{5}{88} & -\frac{270}{1303} & 1 \\ \end{array} \right),\left( \begin{array}{ccccc} 7 & 9 & 6 & 0 & 8 \\ 0 & \frac{20}{7} & -\frac{17}{7} & 10 & -\frac{25}{7} \\ 0 & 0 & -\frac{44}{5} & 34 & -14 \\ 0 & 0 & 0 & -\frac{1303}{88} & \frac{939}{88} \\ 0 & 0 & 0 & 0 & -\frac{4552}{1303} \\ \end{array} \right)$

• "can use Span (;;)" - and thus, Mathematica has caught up to MATLAB. – J. M. is in limbo Oct 24 '15 at 10:54

I believe LUDecomposition is builtin.

LUDecomposition[{{1, 6, 1, 9}, {4, 4, 9, 7}, {10, 4, 10, 2}, {10, 3, 10, 2}}]

{{{1, 6, 1, 9}, {4, -20, 5, -29}, {10, 14/5, -14, -(34/5)}, {10, 57/20, 57/56, 11/7}}, {1, 2, 3, 4}, 1}

• The OP is aware of this, given he compares the times obtained with the builtin. – enzotib Apr 3 '15 at 17:41