# Can anyone give me some suggestions about optimize my Matlab codes?

recently, I try to write a Matlab codes to implement a sparse approximation inverse factorization method proposed by M. Benzi in his paper http://www.mathcs.emory.edu/~benzi/Web_papers/ainv.pdf this codes to generate the incomplete factorization as $A^{-1} = ZD^{-1}W$ or $A^{-1} = \hat{Z}W$. I had written the following codes: then this codes can run with any bugs. However, this codes seems to be very expensive, when I use it to handle a 2961*2961 test matrix, it need one night to product the matrices $W$ and $Z$. So I think that the codes should be optimized to improve the efficiency. So I want to ask someone to give me some suggestions to improve this codes. I also note there are many "for-loop" in my codes, maybe these "for-loop" can be replaced by means of some strategies. Up to now, I have no ideas to improve this codes, so I am turning to everyone.

• In matlab command line, type 'profview' which brings up the Profiler window. Try to locate the bottleneck of your codes! Just a guess by chance, accessing entries of a sparse matrix can be very slow (if it is the case, you can optimise by operating on the [i,j,v] structure of the matrix directly or even use another sparse format. – Hui Zhang Nov 1 '14 at 16:00
• -1 for screenshot instead of copy-and-pasted reusable source. Please change that. – Federico Poloni Nov 1 '14 at 20:07

It is not slow because you've coded the algorithm poorly, or because Matlab is slow at executing it. The algorithm itself is too slow. I'm not quite sure what it does, but I can explain why it's so slow.

In one part of your function you have the following code structure:

for i=1:n
for j=i+1:n
for k=1:i
<body>
end
end
end


This triple nested loop executes its body $$\frac16n(n^2-1)$$ times, so the time complexity of the algorithm is $O(n^3)$. For a matrix of size $n=3000$, this is $\approx 5\times 10^{9}$ times that the body of that loop alone is executed. Thus this algorithm is prohibitively expensive for a matrix of that size.

As for why this algorithm has complexity $O(n^3)$, I don't know. I don't really see where it uses sparsity, because it accesses each element $w_{k,j}$ in that loop, although $w$ starts as speye, so I suspect it's meant to be sparse. Typically an algorithm for a sparse matrix would try to keep the number of operations at something like $O(n)$, not $O(n^3)$. Does the original paper definitely specify an $O(n^3)$ algorithm, or does it have a lower complexity?

If you replace the for-loops with some other looping construct, the number of iterations will stay the same, and this is unlikely to help very much, unless there is nothing else that can be done. The right thing to do is to figure out how to bring the time complexity down (to $O(n^2)$ or lower), but that is a question of algorithm design, and is also a completely different issue from optimizing Matlab code. It is also important to actually take advantage of sparsity, and I don't see how this algorithm does that now.

• Thank you for your detailed analysis. In fact, there exists the FORTRAN-like codes of this method, unfortunately, I cannot read the FORTRAN codes. The numerical examples (in some papers) using this kind of FORTRAN codes show that their algorithm is still very efficient. I can only use the Matlab program, so that is why I turn to the Matlab codes. – Hsien-Ming Ku Nov 1 '14 at 21:48
• @Hsien-MingKu The time complexity of the algorithm is the issue, not whether it's Fortran or Matlab. – Kirill Nov 1 '14 at 22:19