2

I have a CHOLMOD factorization of a sparse matrix H, and I want to edit the sparse representation of the upper, lower, and block diagonal factors. How can I do this? When I run the below, the last line doesn't work.

H = sprand(10,10,0.5)
fac = ldltfact(H; shift=0.0)
fD = fac[:D]
D = Base.SparseArrays.CHOLMOD.Sparse(fD)

And is there any way to go in the reverse direction from a sparse matrix to a CHOLMOD.factor?

jjjjjj
  • 1,152
  • 1
  • 14
  • 30

1 Answers1

2

Extracting the relevant factorization matrices of ldltfact can be a little tedious. The following example shows an example similar to the one in the question with a final test that the extracted matrices recover the original factorized one:

srand(1)
pre = sprand(10,10,0.5)
H = pre + pre' + speye(10,10)

fac = ldltfact(H; shift=0.0)
P = sparse(1:size(H,1),fac[:p],ones(size(H,1)))
LD = sparse(fac[:LD]) # this matrix contains both D and L embedded in it

L = copy(LD)
for i=1:size(L,1)
  L[i,i] = 1.0
end

D = sparse(1:size(L,1),1:size(L,1),diag(LD))

PHP = P*H*P'
LDL = L*D*L'

using Base.Test
@test PHP ≈ LDL

The expected output (and actual on Julia v0.6.3):

julia> @test PHP ≈ LDL
Test Passed

Hope this helps.

Dan Getz
  • 17,002
  • 2
  • 23
  • 41
  • Thanks, this is very helpful! Is there any way to extract the block diagonal factor `D` in sparse form from `LD`? I as because I need to perform a positive definite modification to the 0 and negative eigenvalues, and I don't want to change the whole matrix by performing the factorization with an entire diagonal modification. – jjjjjj Jan 29 '18 at 23:21
  • Suppose that I've edited the sparse representation of `LD`. Is there any nice way to convert it back to a `CHOLMOD.factor`? (per edit) – jjjjjj Jan 30 '18 at 04:07
  • 1
    @jjjjjj Don't think it is easy to convert the edited `LD` into a CHOLMOD.factor, since the functions are basically wrappers around the CHOLMOD library interface. But with the factorization available, the linear algebra should be faster without getting back to use CHOLMOD. The fast `solve` is defined for triangular matrices also. – Dan Getz Jan 30 '18 at 10:31