Operator norm (p = 2) for sparse matrices is not implemented

[goggle]

Using opnorm on a sparse matrix (of type SparseMatrixCSC) gives this error:

ERROR: ArgumentError: 2-norm not yet implemented for sparse matrices. Try opnorm(Array(A)) or opnorm(A, p) where p=1 or Inf.

Here is a piece of code which produces that error:

using LinearAlgebra
using SparseArrays
A = SparseMatrixCSC([1.0 2.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.3])
opnorm(A)

I am interested in working on that issue, but I would need some guidance.

The proper way to solve this (in my opinion) is to implement svdvals for sparse matrices, so that opnorm can be modifed to return the largest singular value of A by using svdvals for sparce matrices. If this is the way to go, which algorithm should be used to implement svdvals for sparse matrices? There have been many discussions about implementing SVD for sparse matrices (e.g. JuliaLang/julia#6610) a long time ago...

If the proposed solution is not the desired way to resolve this, why can't we just replace the error message by returning opnorm(Array(A)) as it is suggested in the error message itself?

[dkarrasch]

Returning opnorm(Array(A)) silently and by default might be really bad because of massive memory consumption, and should be done consciously by the user. I think that's the rationale for the error message.

[goggle]

@dkarrasch I agree, silently returning opnorm(Array(A)) is not a good solution.

But what are Julia's plans to resolve this issue? To resolve this issue in Julia, I assume we would need to have something like a Lanczos Bidiagonlization algorithm implemented. Is this the way to go or rather out-of-scope and should be kept in external Julia packages?

[dkarrasch]

I'm no authority here, but I guess this is a bit out of scope to provide one standard way to do it. You could proceed as you said, compute the topmost singular value. Now here is exactly the point: which method/implementation do you want to use? I'm not sure one can solve that by means of stdlib-functions (which kind of implies that this is not going to be done in SparseArrays). You could use Arpack.jl, or perhaps something from IterativeSolvers.jl, or TruncatedSVD.jl. In fact, you may even appreciate the fact that you can choose a tailored package to solve your specific problem...

[goggle]

Since sparse matrices are part of stdlib, shouldn't they be as feature-complete as dense arrays?

It is currently possible to solve linear systems A * x = b where A is a sparse matrix. I took a quick look at the code, and it seems to use LU factorization to achieve that. Although this is not the preferred way to solve such problems in many real-world applications, the Julia standard library offers a solution that can be good enough for many use cases. This does not mean that it is the only way to solve such a linear system of equations, it's just the way that the Julia standard library provides. So in my opinion, the Julia standard library should provide a feature-full implementation of sparse arrays, which includes an implementation of SVD, opnorm(A, 2), rank, etc. Otherwise one could argue, that there is no need for sparse arrays in stdlib at all, why not use an external package for that?

It would be really interesting to hear the opinions from several Julia maintainers.

[StefanKarpinski]

Otherwise one could argue, that there is no need for sparse arrays in stdlib at all, why not use an external package for that?

I very much want to move sparse arrays out of stdlib.

[goggle]

Interesting.

I guess having SparseArrays in a separate package would make it much easier to add features (like SVD, solving linear systems with an iterative solver (maybe by using something from IterativeSolvers.jl), adding more common matrix decompositions).

On the other hand, sparse matrices are a fundamental data structure which are used in many different applications. Not having them in the standard library would probably worsen the Julia out-of-the-box experience for many users...

[ViralBShah]

There are lots of packages that have fantastic experience for users. I think users will be better served with a faster moving sparse matrix package that lives outside, where capabilities can be added quickly, and new experiments can be carried out.

We'll move Sparse and SuiteSparse out whenever the stdlib versioning stuff is figured out.

[rayegun]

@ViralBShah is this out of scope for SparseArrays even once we move it out?

[ViralBShah]

Not sure I understand. You mean implementing the opnorm? We should implement it all in this package for now.