Error in A_ldiv_B! with sparse matrix

[CodeLenz]

Hi.

I just found a bug in A_ldiv_B! when used with a sparse matrix
factored with lufact!

The error is trigered with

A = spzeros(3,3)
A[1,1] = 1
A[2,2] = 1
A[3,3] = 1
b = ones(3)

lufact!(A)

A\b  --> trigs the error.

ERROR: BoundsError()
in istriu at sparse/sparsematrix.jl:1748
in unsafe_copy! at array.jl:41
in A_ldiv_B! at linalg/sparse.jl:209

julia> versioninfo()
Julia Version 0.3.0-prerelease+3952
Commit b2fd3af (2014-06-30 01:39 UTC)
Platform Info:
System: Windows (x86_64-w64-mingw32)
CPU: Intel(R) Core(TM) i5-3320M CPU @ 2.60GHz
WORD_SIZE: 64
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)
LAPACK: libopenblas
LIBM: libopenlibm

It was working with the previous versions.

Sincerely yours,
Eduardo.

[ViralBShah]

@JeffBezanson It was strange that the above commit fixed this bug. I wonder if this is a compiler issue.

[tkelman]

And this wasn't caught by any of the existing tests in test/sparse.jl, at least not on my machine (and I can reproduce the error at the same commit). Wouldn't hurt to throw in another, right?

[andreasnoack]

@ViralBShah I don't think that the commit changes anything for this example. With latest master I still get the BoundError.

@CodeLenz I am not sure it is a bug. When calling lufact! you destroy the input matrix which is signalled by the exclamation mark. The row and column pointers are transformed to zero base and therefore not valid in Julia anymore.

What you might want to do is either

A = speye(3)
b = ones(3)
Af = lufact!(A)
Af\b

or

A = speye(3)
b = ones(3)
lufact(A) # No exclamation mark
A\b

[CodeLenz]

Hi. It worked with the previous versions. As I understand, lufact!(A) stores the LU factors of
A in place and A\b makes use of it to solve the system. But I am new to Julia, so if its my mistake, I am realy sorry for the noise.

[jiahao]

lufact!(A) stores the LU factors of A in place and A\b makes use of it to solve the system.

If you mean to factorize A and use the factorization to solve a linear problem, you'll need to save the factorization object and call backslash on that.

julia> B=copy(A); Z=lufact!(B)
UMFPACK LU Factorization of a 3-by-3 sparse matrix
Ptr{Void} @0x00007f85ab78d8c0

julia> typeof(B) #Julia still thinks it is a sparse matrix
SparseMatrixCSC{Float64,Int64} (constructor with 2 methods)

julia> Z\b
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> @which Z\b #Uses the factorization object
\{TF<:Number,TB<:Number,N}(F::Factorization{TF<:Number},B::AbstractArray{TB<:Number,N}) at linalg/factorization.jl:798

julia> @which A\b #Doesn't use the factorization object
@\{TA,TB,N}(A::AbstractArray{TA,2},B::AbstractArray{TB,N}) at linalg/generic.jl:227

A\b will solve the system without using the factorization. But since lufact!(A) destroys A, there is no guarantee that A represents a valid Julia object, and so calling A\b can fail.

[CodeLenz]

@jiahao I see ...thanks. I will look more carefully before posting another bug report. But I have no idea why it used to work before and the ! mark lead me to this conclusion....so....my bad.

Thank you all for this amazing tool. I just found it and I am porting all my tools (FInite element and optimization) to it.

[jiahao]

But I have no idea why it used to work before

It could be that you were just lucky before, since the specific memory access patterns can vary a lot.

[JeffBezanson]

The documentation could use some extra clarification here. In the vast majority of cases, f!(A) means that the result is written into A. Instead, this is a strange kind of function that returns a result but destroys its input in the process.

[jiahao]

Is there a way to make it an error to access A after calling lufact!(A)? I honestly have no idea what UMFPACK does to A.

[ViralBShah]

I do think that we should not allow lufact!(A) where the trashed input is then unusable. @dmbates What do you think? This seems like it will cause more users to trip up, expecting dense matrix like behaviour. I don't think that documentation is sufficient for this pitfall.

[andreasnoack]

All of the xfact! functions work something like this. For the dense methods, the input matrix stores information about the factorisation, but it is only useful if used through the factorisation object, not the input matrix.

[ViralBShah]

Yes, I guess that even though the solution is returned in the input, it is not really in usable form.

[mlubin]

How about a test for valid indices when calling A\b?

[dmbates]

A simple test for valid indices is whether the first element of the colptr member is 0 or 1. To be a valid Julia matrix it should be 1. To be a valid SuiteSparse matrix it should be 0.

It turns out that the only difference in behavior between lufact and lufact! with regard to the sparse matrix is whether the function decrements the indices of the original matrix in place. I think I would agree with @ViralBShah that the best thing to do is to remove the lufact! method for sparse matrices. The sparse case is different from the dense case in that you can't overwrite a sparse matrix with its factorization in most cases. Thus the lufact! method doesn't make sense unless you count saving the allocations of the index vectors. I initially thought that would be helpful if the sparse matrix was the result of other functions and was not going to be named but, as we have seen, it just leads to confusion.

[tkelman]

The right (medium-level, with reusable factorizations - A\b works fine now as the high level way) way to eventually do sparse is separate symbolic and numerical factorizations. symfact only looks at sparsity structure and returns a preallocated object with enough predicted space needed for numfact!. Something like this isn't exposed yet in the umfpack interfaces, is it?

[ViralBShah]

I am pretty sure such a thing will be there in the umfpack interface.

[tkelman]

Whoops, yep, there they are in umfpack.jl. Maybe an area for future refactoring and generalization to expose those, I think it could be doable to present a unified medium-level API that all of umfpack, cholmod, pardiso, mumps, hsl, etc could work through. I think pardiso and mumps even let you specify whether your indices are 0-based or 1-based, so the input formatting question is different.

It could be almost-as-dangerous (maybe more since it doesn't error?) and confusing to accidentally misuse a dense matrix that has had lufact! called on it. For dense it's not as obvious that the array now contains factorization data since there isn't the same validity check we get with 0-vs-1 indexing in sparse. lufact! for sparse does save a little bit of copying if you don't plan on using the input matrix again while making calls to umfpack, plus it's an easily reversible operation to turn it back into a valid 1-based Julia sparse matrix (assuming you're aware of exactly what's going on behind the scenes).

As @dmbates said, "in place" is not so much a thing with sparse. You do often want multiple factorizations in a row with the same sparsity pattern but changed nzval.

[ViralBShah]

Yeah, I think this is the right thing to do. It is certainly doable in 0.3.

[Jutho]

All the xfact! methods potentially destroy the input. It is good to have them when you can afford to loose the original matrix, sparse matrix, but like a lot of the ! functions, even those that do not destroy any input (transpose!, permutedims!, A_mul_B! family), it might be possible not to export them, so as not to confuse new users.

[mlubin]

I think cheap error checking is an easy way to avoid confusion. Removing the ! methods also breaks generic code that can work with both dense and sparse matrices.

[jiahao]

Another possibility is to simply nudge the pointer values back to 1-based indexing after calling UMFPACK functions.

[dmbates]

@jiahao Seems reasonable.

[ViralBShah]

We have now removed ! methods for sparse factorizations.

[ViralBShah]

@andreasnoack Should we disallow cholfact! too? It gives a DomainError currently:

julia> cholfact!(sparse(rand(5,5)))


ERROR: DomainError:

 in chol! at linalg/cholesky.jl:133
 in chol! at linalg/cholesky.jl:71
 in cholfact! at linalg/cholesky.jl:98

[ViralBShah]

I see that it has a fallback to the generic implementation. Perhaps we need an explicit method to disallow it for sparse?

[andreasnoack]

I think it is safe to restrict that method to StridedMatrix. That should fix the issue, so I'll do that.

[andreasnoack]

Fixed in 9c72dec

[tkelman]

Removing the ! methods also breaks generic code that can work with both dense and sparse matrices.

This is still a potential concern. But I don't think our current cholfact! API is really generic enough to cover sparse effectively in its current form. What would probably be better for dense-sparse generality is deciding on (then implementing, testing, and documenting it properly) a good medium-level API for reusable sparse factorizations, something like symfact and numfact!. Dense can be expressed as the special-case performance optimization here, within a framework that's generic enough to do sparse correctly. symfact for dense just happens to be a no-op (or identity, returning the same array space).

[andreasnoack]

I agree. The generic cholfact! computes the factorization correctly, but is not useful for sparse matrices. It could be interesting to experiment with sparse factorizations in Julia, but we should probably do the experiments outside base.

[ViralBShah]

In case you did not already know, see: https://github.com/JuliaSparse/MultiFrontalCholesky.jl

[tkelman]

That's great to see a start for a pure-Julia sparse solver implementation. So far that's working at a pretty low level though, the API refinement will need to happen somewhat separately from and in parallel with the pure-Julia solver implementations. Ideally we want a design somewhat like MathProgBase where different low-level solvers (some in Julia, some in C or C++ or Fortran) can plug into a uniform set of interfaces, to present the same API to consumers (whether users directly working with sparse linear algebra, or higher-level libraries like optimization solvers).

This probably should happen outside of base to start with. Though to truly resolve JuliaLang/LinearAlgebra.jl#136 and make sparse matrices first-class citizens, the way base Julia deals with dense matrices will eventually have to adapt in order to generically handle either dense or sparse matrices for all operations.