Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am trying to update some parts of an specific matrix as rapidly as possible. In what follows, I first set up the basics things that I want to use

d = 300;
A = RandomReal[20, {d, d}];
n = Dimensions[A][[1]];
i = 1; num = 2;
Id = SparseArray[{{k_, k_} -> 1.}, {n, n}];
LU = LinearSolve[A];
mat = (1./(Norm[A, 2]^2))*ConjugateTranspose[A];
mat // MatrixPlot

wherein Id, denotes the Identity matrix, A is a dense input matrix, and then I wish to update the matrix mat, (which is an approximate inverse of the matrix A), by entering num=2, columns of the exact inverse that will be obtained by solving two linear systems (using LinearSolveFunction). To this end, I use the following piece of code, in which after obtaining the num=2 columns of the matrix inverse, they must be replaced as the first and the second columns of the matrix mat at the end of each cycle of While: (please forgive me, if I write the codes in a very rough way!)

While[i <= num,
  {ll = Id[[All, i]];
   ith = Chop@LU[ll];
   mat[[All, i]] = ith;
mat // MatrixPlot

Considering the above dense matrix A, it works and can update the columns of the matrix mat. My problem is here, if I use a sparse matrix, then for low dimensions it works rapidly, while for higher dimensions it takes too much time to update the columns of the matrix mat. I mean, if we use the following matrix

A = SparseArray[{{i_, i_} -> 
 RandomReal[3], {i_, j_} /; Abs[i - j] == 1 -> RandomReal[2],
{i_, j_} /; Abs[i - j] == 8 -> RandomReal[1]}, {d, d}];

while d=3000. I would like to ask you experts about that: how we could accelerate (in terms of computational time) this process for dense and sparse matrices in a simple uniform piece of code?

Any tip or help will be cheerfully thanked.

share|improve this question
Useful references: (1) (2) (3) -- you'll notice these are all from Leonid Shifrin. – Mr.Wizard Feb 24 '13 at 13:59
I would add this link, where I outlined the general efficiency problem with by-element updates of sparse matrices. You would need a sparse array implementation based on something like binary trees to make updates efficient, I think. – Leonid Shifrin Feb 24 '13 at 15:18
I think the above piece of codes fail even when we test a 300*300 sparse matrix. I mean sometimes in my MMA 8, it gives the results and sometimes when I increase $num$ to 3 (for instance), it fails and the MMA generates a beep! – Fazlollah Soleymani Feb 24 '13 at 16:22
In my idea, the best way is to extract the diagonal elements of the matrix inverse, and then update your matrix $mat$. In such a way, the norm of the matrix $Id-A.mat$ will decrease much more and with one replacement you might obtain the best possible approximate inverse. The only problem is that how to extract or find the diagonal entries of the matrix inverse very fast for a very large sparse matrix! – Fazlollah Soleymani Feb 24 '13 at 16:26

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.