factorize is not implemented for SharedSparseMatrix

[ingenieroariel]

Using Julia's master (from Feb 17, 2014) on OSX:

addprocs(6)
using ParallelSparseMatMul
m,n,p = 100,30,.1
S = shsprandn(m,n,p)
S = factorize(S)
ERROR: no method factorize(SharedSparseMatrixCSC{Float64,Int64})

[ingenieroariel]

Tried the following workaround without success:

julia> cholfact(localize(shsprandn(100,30,.1)));
CHOLMOD error: argument missing
ERROR: CholmodException()

As a comparison:

julia> cholfact(sparse(randn(100,30)))
CHOLMOD warning: matrix not positive definite

CHOLMOD factor:  :  100-by-100 not positive definite (column 30)
  scalar types: SuiteSparse_long, real, double
  supernodal, LL'.
  ordering method used: AMD
         0:0
         1:1
         2:2
         3:3
         4:4
         5:5
         6:6
         7:7
    ...
        96:96
        97:97
        98:98
        99:99
  col: 0 colcount: 100
  col: 1 colcount: 99
  col: 2 colcount: 98
  col: 3 colcount: 97
  col: 4 colcount: 96
  col: 5 colcount: 95
  col: 6 colcount: 94
  col: 7 colcount: 93
    ...
  col: 96 colcount: 4
  col: 97 colcount: 3
  col: 98 colcount: 2
  col: 99 colcount: 1
  ssize 100 xsize 10000 maxcsize 1 maxesize 1
  supernode 0, col 0 to 99. nz in first col: 100.
  values start 0, end 10000
  col 0
         0: 6.0464
         1: 1.5983
         2: -0.58295
         3: 1.2476
    ...
  col 97
        97: 0
        98: 0
        99: 0
  col 98
        98: 0
        99: 0
  col 99
        99: 0
  nz 5050  OK

[madeleineudell]

B = sprandn(100,30,.1); cholfact(localize(share(B))) works. I believe the reason that your workaround above doesn't work is that shsprandn doesn't guarantee that entries have unique values; ie I can have an entry [3,4] = .4 and [3,5] = -.2. This is fine if I'm only multiplying by the matrix (the entries are just added) but makes no sense for factorizations. I'll open a new issue for it.

[ingenieroariel]

The original ticket was about overriding factorize in ParallelSparseMatMul to support SharedSparseMatrix. Would a parallel factorization make any sense or is it one of the operations that is better done serialized?

[ingenieroariel]

For future readers, I think @madeleineudell meant:

B = sprandn(100, 30, 0.1); cholfact(localize(share(B)))

and the new issue she created is:

https://github.com/madeleineudell/ParallelSparseMatMul.jl/issues/3

[madeleineudell]

cholfact call suitesparse, which already has a parallel factorization in it. i'm not sure if julia's cholfact calls it by default or not. but parallel cholesky is a hard enough problem that it's not a great idea to reimplement it: see http://beowulf.csail.mit.edu/18.337-2012/projects/ParallelLinearAlgebrainJulia.pdf for a noble attempt that failed. ParallelSparseMatMul is really a package intended to support multiplication by A and A'. You can solve A\b using eg lsqr https://github.com/JuliaLang/IterativeSolvers.jl by calling only multiplication by A and A', and that's the direction this project is going. It's an "operator" approach to least squares, rather than an element-conscious factorization approach.

On Tue, Feb 18, 2014 at 10:00 AM, Ariel Núñez notifications@github.comwrote:

The original ticket was about overriding factorize in ParallelSparseMatMul to support SharedSparseMatrix. Would a parallel factorization make any sense or is it one of the operations that is better done serialized?

Reply to this email directly or view it on GitHubhttps://github.com/madeleineudell/ParallelSparseMatMul.jl/issues/2#issuecomment-35413372 .

Madeleine Udell PhD Candidate in Computational and Mathematical Engineering Stanford University www.stanford.edu/~udell

[ingenieroariel]

Thanks a lot @madeleineudell ! I will look into IterativeSolvers.jl and try to write a sample program in Julia that solves a simple Lasso.

I am getting ahead of myself by looking into ParallelSparseMatMul before learning Julia well enough (I code in Python) but your post to the mailing list about the performance made me want to play with it right away.

[madeleineudell]

Ariel, you can also try this package (expect lots of bugs, though): https://github.com/madeleineudell/ParallelSparseRegression.jl

On Tue, Feb 18, 2014 at 10:15 AM, Ariel Núñez notifications@github.comwrote:

Thanks a lot @madeleineudell https://github.com/madeleineudell ! I will look into IterativeSolvers.jl and try to write a sample program in Julia that solves a simple Lasso.

I am getting ahead of myself by looking into ParallelSparseMatMul before learning Julia well enough (I code in Python) but your post to the mailing list about the performance made me want to play with it right away.

Reply to this email directly or view it on GitHubhttps://github.com/madeleineudell/ParallelSparseMatMul.jl/issues/2#issuecomment-35414980 .

Madeleine Udell PhD Candidate in Computational and Mathematical Engineering Stanford University www.stanford.edu/~udell