Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Solving with sparse cholesky returns two-dimensional array #28985

Closed
clason opened this issue Aug 31, 2018 · 10 comments
Closed

Solving with sparse cholesky returns two-dimensional array #28985

clason opened this issue Aug 31, 2018 · 10 comments
Labels
domain:arrays:sparse Sparse arrays domain:linear algebra Linear algebra kind:bug Indicates an unexpected problem or unintended behavior

Comments

@clason
Copy link

clason commented Aug 31, 2018

In v1.0, using the Cholesky decomposition of a sparse array to solve a linear system for the transpose results in a two-dimensional array instead of a vector:

using LinearAlgebra,SparseArrays
N = 4
A = sparse(I,N,N)
C = cholesky(A)
b = randn(N)
C'\b

returns

4×1 Array{Float64,2}:
...

By contrast, C\b (correctly) returns

4-element Array{Float64,1}:
...

(The corresponding test for Array(I,N,N) throws a ERROR: MethodError: no method matching adjoint(::Cholesky{Float64,Array{Float64,2}}).)

@ViralBShah ViralBShah added the domain:arrays:sparse Sparse arrays label Aug 31, 2018
@StefanKarpinski
Copy link
Sponsor Member

@Sacha0, @andreasnoack?

@andreasnoack andreasnoack added the kind:bug Indicates an unexpected problem or unintended behavior label Sep 10, 2018
@andreasnoack
Copy link
Member

This is a bug and should be fixed.

@ViralBShah
Copy link
Member

Yes, we must fix this. The culprit is probably in here:

\(adjL::Adjoint{<:Any,<:Factor}, B::Dense) = (L = adjL.parent; solve(CHOLMOD_A, L, B))
\(adjL::Adjoint{<:Any,<:Factor}, B::VecOrMat) = (L = adjL.parent; Matrix(solve(CHOLMOD_A, L, Dense(B))))

@raghav9-97
Copy link
Contributor

I would like to work on this if nobody is currently working on it.

@ViralBShah
Copy link
Member

Go for it.

@raghav9-97
Copy link
Contributor

I think the problem is explicit conversion to Matrix in \(adjL::Adjoint{<:Any,<:Factor}, B::VecOrMat) = (L = adjL.parent; Matrix(solve(CHOLMOD_A, L, Dense(B))))
Is that correct?

@ViralBShah
Copy link
Member

I believe the matrix and the vector cases need to be separated out.

@raghav9-97
Copy link
Contributor

raghav9-97 commented Dec 17, 2018

\(adjL::Adjoint{<:Any,<:Factor}, B::Vector) = (L = adjL.parent; Vector(solve(CHOLMOD_A, L, Dense(B))))
\(adjL::Adjoint{<:Any,<:Factor}, B::Matrix) = (L = adjL.parent; Matrix(solve(CHOLMOD_A, L, Dense(B))))

Something like this??

@ViralBShah
Copy link
Member

Probably. It needs to be tried out and all tests should pass.

@andreasnoack
Copy link
Member

Something like this??

Yes that should work. While you are at it, please allow that B can be StridedVector/StridedMatrix. Also, it's slightly nicer to use lower case for vectors (but still upper case for matrices). Finally, please also split the function over multiple lines instead of using ;.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
domain:arrays:sparse Sparse arrays domain:linear algebra Linear algebra kind:bug Indicates an unexpected problem or unintended behavior
Projects
None yet
Development

No branches or pull requests

5 participants