COO to CSC format with the same I,J vectors

Question

Right now, I am calling K = sparse(I,J,V,n,n) function to create sparse (symmetric) K matrix in Julia. And, I am doing this for many steps.

Due to memory and efficiency considerations, I would like to modify the K.nzval values, instead of creating a new sparse K matrix. Note that the I and J vectors are the same for each step, but the non-zero values (V) are changing at each step. Basically, we can say we know the sparsity pattern in COO format. (I and J are not ordered and may have multiple (I[i],J[i]) entries)

I tried to order my COO format vectors to relate to the CSC/CSR format storage. However, I found it non-trivial (at least, for now).

Is there a way to do this or a magical "sparse!" function? Thanks,

Here is an example code that is relevant to my question.

n=19 # this is much bigger in reality ~ 100000. It is the dimension of a global stiffness matrix in finite element method, and it is highly sparse! 
I = rand(1:n,12)
J = rand(1:n,12)
#
for k=365
  I,J,val = computeVal()  # I,J are the same as before, val is different, and might have duplicates in it. 
  K = sparse(I,J,val,19,19)
  # compute eigs(K,...)
end

# instead I would like to decrease the memory/cost of these operations with following
# we know I,J
for k=365
  I,J,val = computeVal()  # I,J are the same as before, val is different, and might have duplicates in it. 
  # note that nonzeros(K) and val might have different size due to dublicate entries.
  magical_sparse!(K,val)
  # compute eigs(K,...)
end

# what I want to implement 
function magical_sparse!(K::SparseMatrixCSC,val::Vector{Float64}) #(Note that this is not a complete function)
    # modify K 
    K.nzval[some_array] = val
end

edit:

A more specific example is given here.

n=4 # dimension of sparse K matrix
I = [1,1,2,2,3,3,4,4,1,4,1]
J = [1,2,1,2,3,4,4,3,2,4,2]
# note that the (I,J) -> (1,2) and (4,4) are duplicates.
V = [1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.]

function computeVal!(V)
  # dummy function
  # return modified V
  rand!(V) # this part is involed, so I will just use rand to represent that we compute new values at each step for V vector. 
end

for k=1:365
  computeVal!(V)
  K = sparse(I,J,V,n,n)
  # do things with K
end

# Things to notice:
# println(length(V))       -> 11
# println(length(K.nzval)) -> 8

# I don't want to call sparse function at each step. 
# instead I would like to decrease the cost of these operations with following
# we know I,J
for k=1:365
  computeVal!(V)
  magical_sparse!(K,V)
  # do things with K
end

# what I want to implement 
function magical_sparse!(K::SparseMatrixCSC,V::Vector{Float64}) #(Note that this is not a complete function)
  # modify nonzeros of K and return K  
end

It would help to have more details about the COO format you have (ideally with some pasted code as well) — Dan Getz, Dec 19 '17 at 21:21
Use `nonzeros(K) .= val` to implement the magical function or even directly randomize `nonzeros(K)` using `rand!` — Dan Getz, Dec 20 '17 at 05:35
That function helps! Thank you. But, there are still parts that need to be worked on. I also edited the question to clarify some lines. In fact, `val` is not a random vector. It is the array corresponding to new computed values for `[I,J]` coordinates. I guess one needs to sum over the duplicates to be able to use your suggestion! I will give it a try! Thanks again. — berkin, Dec 20 '17 at 19:21
for example; `[I,J,V]` might have `[1,1,0.1]` and `[1,1,0.6]` as entries. In this case, `sparse` function sums them over the same `I`,`J` coordinate, eg. `[1,1,0.7]`. — berkin, Dec 20 '17 at 20:37
If `I`, `J`, `V` are vectors, then `[1,1,0.1]` is not an entry/element. What is an example of the output of `computeVal()`? — Dan Getz, Dec 20 '17 at 21:19
Will every entry in the sparse matrix appear in `I,J` or could some entries of the sparse matrix be missing from the update? If so, do they need to be zeroed or left as before? — Dan Getz, Dec 20 '17 at 21:39
Yes, `I`,`J`,`V` are vectors. `[I,J,V]` is a matrix with `I`,`J`,`V` as columns. `computeVal()` returns `I`,`J`,`V` separately (as vectors), but I and J are the same for all the steps (you can argue that the function `computeVal()` does not have to return `I` and `J`, because they are not modified. The function only computes new `V` in each step). Entries in `I` and `J` vectors correspond to row and column index of the corresponding nonzero entry, respectively. All the entries in sparse matrix will appear in `I` and `J` at every step as nonzeros. — berkin, Dec 20 '17 at 21:53
Before optimizing the solution, does `foldl((x,y)->(x[y[1],y[2]]+=y[3];x),fill!(nonzeros(K),0.0), zip(I,J,val))` work? It may a bit cryptic, but try it after setting `I,J,val`. — Dan Getz, Dec 20 '17 at 22:16
Got dimension errors! (But I added a specific example where you can run the algorithm yourself, see the edit). Thanks again! — berkin, Dec 20 '17 at 23:16
Updated the answer for the new setup in the question. It should work, and quite efficiently, but my money is on the question changing a bit more in response ;) — Dan Getz, Dec 21 '17 at 16:28
That's what the magical sparse function is! Thank you very much! Thanks for the patient about the edits on the question :) — berkin, Dec 21 '17 at 21:41
BTW because the indices accessed in `K` are constant between rounds, there is obvious room for optimization by not recalculating the position of changes in `nonzeros(K)`. But if later processing includes eigenvalue calculation, it may be an overkill to consider this now. — Dan Getz, Dec 23 '17 at 08:36

Dan Getz · Accepted Answer · 2017-12-21T17:29:08.927

UPDATE FOR CURRENT QUESTION

According to the changes in the question, the new solution is:

for k=365
  computeVal!(V)
  foldl((x,y)->(x[y[1],y[2]]+=y[3];x),fill!(K, 0.0), zip(I,J,V))
  # do things with K
end

This solution uses some trickery, such as, zeroing K using fill!, which by default returns K which is then used as an initial value for the foldl. Again, using ?foldl should clear up what is going on here.

ANSWER TO OLD QUESTION

Replacing

for k=365
  val = rand(12)
  magical_sparse!(K,val)
  # compute eigs(K,...)
end

with

for k=365
  rand!(nonzeros(K))
  # compute eigs(K,...)
end

should do the trick.

Use ?rand! and ?nonzeros for help on the respective functions.

COO to CSC format with the same I,J vectors

1 Answers1