# How to accelerate updating some parts of sparse matrices?

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

Clear["Global`*"]
SeedRandom[1234];
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;
i++}
];
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.

• 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 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 Feb 24 '13 at 16:26

(*Ksp is the sparse matrix, add a_List to the diagonal of Ksp*)
Ksp += SparseArray[{{i_, i_} :> a[[i]]}, Length[a]];

Two methods to add a list to the diagonal of a sparse matrix are given as below. The performance varies greatly.

n = 100000;

a = RandomReal[1, n];

Ksp = SparseArray[{{i_, i_} :> a[[i]]}, Length[a]];

Ksp2 = Ksp;

(Method 1)

In[42]:= t1 = Timing[Do[Ksp[[i, i]] += a[[i]], {i, Length[a]}]][[1]]

Out[42]= 33.7742

(Method 2)

In[43]:= t2 = Timing[Ksp2 += SparseArray[{{i_, i_} :> a[[i]]}, Length[a]]][[1]]

Out[43]= 0.0312002

In[44]:= t1/t2

Out[44]= 1082.5

• This is not a self-contained answer: what are the values of Ksp and a? I cannot find these in the original post. And why not use DiagonalMatrix[a] instead? – István Zachar Nov 5 '18 at 17:38
• Agreed. Please fill in the details. – Henrik Schumacher Nov 5 '18 at 18:31
• You can't use DiagonalMatrix[a] when Length[a] is large. DiagonalMatrix will create a dense matrix rather than a sparse matrix. – HL REN Nov 5 '18 at 19:59
• Actually, you can: DiagonalMatrix[a, 0, Length[a], SparseArray] (for you too, @István) – J. M. will be back soon Jan 5 at 4:51