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

Rewrite spdiagm with pairs #23757

Merged
merged 1 commit into from
Oct 11, 2017
Merged

Rewrite spdiagm with pairs #23757

merged 1 commit into from
Oct 11, 2017

Conversation

fredrikekre
Copy link
Member

@fredrikekre fredrikekre commented Sep 18, 2017

I think this is nicer than the previous spdiagm(tuple_of_vectors, tuple_of_diags) API. And the same should also be implemented for diagm. See JuliaLang/LinearAlgebra.jl#457 and specifically JuliaLang/LinearAlgebra.jl#457

This also deals with JuliaLang/LinearAlgebra.jl#457 such that we always return a square matrix.

Some things to decide:

  • Should the pair be vector => diagonal or diagonal => vector? I think vector => diagonal is nicer.
  • Should we have a convenience constructor for the single vector case? spdiagm(v::AbstractVector) = spdiagm(0 => v). I imagine that case being quite common and it might be nice to have that?

Still need to implement proper deprecations and such.

@fredrikekre
Copy link
Member Author

Actually, the single vector case could better be written as sparse(Diagonal(v)), that should also be faster.

@fredrikekre
Copy link
Member Author

Here are some datapoints that support the changes here, specifically the removal of the option to supply m, n as the final size. In all registered packages there are ≈120 uses of spdiagm as follows:

  • 69 % are simply spdiagm(x::AbstractVector)
  • 21 % are calls to spdiagm(x::AbstractVector, k, m, n) where m, n are needed to obtain a square matrix.
    julia> spdiagm(rand(2), -1)
    3×2 SparseMatrixCSC{Float64,Int64} with 2 stored entries:
      [2, 1]  =  0.0960447
      [3, 2]  =  0.190498
    
    julia> spdiagm(rand(2), -1, 3, 3) # 3, 3 needed to obtain a square matrix
    3×3 SparseMatrixCSC{Float64,Int64} with 2 stored entries:
      [2, 1]  =  0.64661
      [3, 2]  =  0.119494
    
  • 5.9% are calls with multiple diagonals that will be unaffected of the removal of m, n arguments
  • 4.2 % are calls where m != n, e.g. requesting a non-square matrix.

@deprecate spdiagm(x::AbstractVector) sparse(Diagonal(x))
@deprecate spdiagm(x::AbstractVector, d::Number) spdiagm(x => d)
@deprecate spdiagm(x, d) spdiagm((x[i] => d[i] for i in 1:length(x))...)
@deprecate spdiagm(x, d, m::Integer, n::Integer) spdiagm((x[i] => d[i] for i in 1:length(x))...)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IIUC, the latter now always yields a square matrix whereas the former yields an m-by-n matrix?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yea. I should rewrite these deprecations such that this change won't be breaking.

i += numel
end

return (I,J,V)
return I, J, V
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The former descriptive variable names made this code a bit easier to follow :).

one diagonal, `B` can be a vector (instead of a tuple) and `d` can be the diagonal position
(instead of a tuple), defaulting to 0 (diagonal). Optionally, `m` and `n` specify the size
of the resulting sparse matrix.
Construct a sparse diagonal matrix from pairs of vectors and diagonals.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"pairs" -> "Pairs"?

@ViralBShah
Copy link
Member

ViralBShah commented Oct 1, 2017

Off topic, I also wish there was a better way to name this and similar functions. It is really annoying to have the sp prefixes.

@StefanKarpinski
Copy link
Member

It is really annoying to have the sp prefixes.

The obvious solution would be something like diagm(SparseMatrixCSC, ...).

@ViralBShah
Copy link
Member

diagm(SparseMatrixCSC, ...) is what I thought of too, but it feels too verbose. One other possibility is using the Sparse module and calling Sparse.diagm.

@fredrikekre
Copy link
Member Author

That's what I proposed in JuliaLang/LinearAlgebra.jl#457 and also what we should probably do, see #11557 and #16740. There is currently no API to extend for user arrays for this kind of thing, e.g. zeros, ones, rand etc. But this is unrelated to the current PR.

@fredrikekre fredrikekre changed the title [WIP] Rewrite spdiagm with pairs Rewrite spdiagm with pairs Oct 3, 2017
@deprecate spdiagm(x::AbstractVector) sparse(Diagonal(x))
function spdiagm(x::AbstractVector, d::Number)
depwarn(string("spdiagm(x::AbstractVector, d::Number) is deprecated, use ",
"spdiagm(x => d) instead, which now return a square matrix. To preserve the old ",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"return" -> "returns"?

end
function spdiagm(x, d)
depwarn(string("spdiagm((x1, x2, ...), (d1, d2, ...)) is deprecated, use ",
"spdiagm(x1 => d1, x2 => d2, ...) instead, which now return a square matrix. ",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"return" => "returns"?

end
function spdiagm(x, d, m::Integer, n::Integer)
depwarn(string("spdiagm((x1, x2, ...), (d1, d2, ...), m, n) is deprecated, use ",
"spdiagm(x1 => d1, x2 => d2, ...) instead, which now return a square matrix. ",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"return" -> "returns"?

Copy link
Member

@Sacha0 Sacha0 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These changes strike me as a step forward, independent of where diagm and friends go eventually. Much thanks for thinking through the design here and in upstream threads @fredrikekre! :)

A minor, subjective point: I find k => v more natural than v => k, particularly with the pairs consistently identified as kv. But I'm sure this impression varies, and practically the difference is likely negligible :).

@fredrikekre
Copy link
Member Author

A minor, subjective point: I find k => v more natural than v => k, particularly with the pairs consistently identified as kv. But I'm sure this impression varies, and practically the difference is likely negligible :).

I read v => k as "assign vector v to diagonal k" which I feel is more natural compared to "assign diagonal k to vector v" for k => v.
The reason it is currently called kv is that I implemented as k => v first and then changed. We should definitely change this to match whatever we decide.

@Sacha0
Copy link
Member

Sacha0 commented Oct 7, 2017

Either orientation sounds good on this end :).

@fredrikekre
Copy link
Member Author

Either orientation sounds good on this end :).

I don't particularly care either. Perhaps other people have some opinions?

@mschauer
Copy link
Contributor

mschauer commented Oct 8, 2017

A matrix can be viewed as an association mapping a number (the k in kth diagonal) to a vector (the kth diagonal). In this light, diagm(k =>d) would be more natural.

@martinholters
Copy link
Member

Also, the same vector might be used for different diagonals, but not vice versa.

instead of using one tuple with diagonals and one tuple with vectors
@fredrikekre
Copy link
Member Author

Should be ready now. Wanna take another look @Sacha0 ?

Copy link
Member

@Sacha0 Sacha0 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lgtm from a rapid skim of this iteration! :)

@KristofferC KristofferC merged commit 5295feb into master Oct 11, 2017
@KristofferC KristofferC deleted the fe/spdiagm-pair branch October 11, 2017 07:35
@stevengj
Copy link
Member

Do we still have to call SparseArrays.spdiagm_internal in order to get a non-square matrix?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
sparse Sparse arrays
Projects
None yet
Development

Successfully merging this pull request may close these issues.

9 participants