4
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ Useful references: (1) (2) (3) -- you'll notice these are all from Leonid Shifrin. $\endgroup$ Commented Feb 24, 2013 at 13:59
  • $\begingroup$ 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. $\endgroup$ Commented Feb 24, 2013 at 15:18
  • 1
    $\begingroup$ 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! $\endgroup$ Commented Feb 24, 2013 at 16:22
  • 1
    $\begingroup$ 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! $\endgroup$ Commented Feb 24, 2013 at 16:26

1 Answer 1

Reset to default
-1
$\begingroup$ 
(*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

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ 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? $\endgroup$ Commented Nov 5, 2018 at 17:38
  • $\begingroup$ Agreed. Please fill in the details. $\endgroup$ Commented Nov 5, 2018 at 18:31
  • $\begingroup$ You can't use DiagonalMatrix[a] when Length[a] is large. DiagonalMatrix will create a dense matrix rather than a sparse matrix. $\endgroup$ Commented Nov 5, 2018 at 19:59
  • 2
    $\begingroup$ Actually, you can: DiagonalMatrix[a, 0, Length[a], SparseArray] (for you too, @István) $\endgroup$ Commented Jan 5, 2019 at 4:51
  • $\begingroup$ Good idea. I agree with you. $\endgroup$ Commented Dec 25, 2019 at 11:23

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.