4

Trying to replicate a calculation from Matlab in Julia but am having trouble converting a single column complex array into a sparse diagonalized array for matrix multiplication.

Here is the Matlab code I am trying to replicate in Julia:

x*diag(sparse(y))

where x is of size 60,600000 and is of type: double, and y is size 600000,1 and is of type: complex double.

phipsgabler
  • 20,535
  • 4
  • 40
  • 60
user4434674
  • 193
  • 1
  • 1
  • 8

2 Answers2

3

You can use Diagonal for that:

using LinearAlgebra
x=rand(6,60)
y=rand(Complex{Float64},60,1)
x*Diagonal(view(y,:))

I have used view(y,:) to convert y to a Vector - this is here a dimension drop operator you can also use shorter form vec(y) instead. Depending on what you want to do you might explicitly say that you want first column by view(y,:,1).

Note that Diagonal is just a sparse representation of a matrix.

julia> Diagonal(1:4)
4×4 Diagonal{Int64,UnitRange{Int64}}:
 1  ⋅  ⋅  ⋅
 ⋅  2  ⋅  ⋅
 ⋅  ⋅  3  ⋅
 ⋅  ⋅  ⋅  4

Another option that might cover more use case scenarios are BandedMatrices:

using BandedMatrices
x*BandedMatrix(0=>view(y,:))

Note that BandedMatrix uses set of pairs for bands, where band 0 is actually the diagonal.

Przemyslaw Szufel
  • 40,002
  • 3
  • 32
  • 62
  • In this particular case, `x .* y'` should result in the same, right?. It would be interesting to benchmark what is better. (And congrats on 10k!) – phipsgabler Aug 24 '20 at 06:31
  • 2
    You can write `Diagonal(vec(y))` instead of `Diagonal(view(y, :))`. `vec` is also non-allocating, but you end up with a nicer type signature, and it's easier to write. `vec` is a really nice function. – DNF Aug 24 '20 at 07:07
  • @phipsgabler Thank you! `x .* y'` seem to return a different result. – Przemyslaw Szufel Aug 24 '20 at 08:52
  • @DNF yeah I should type `vec(y)` since it is shorter but somehow every time I automatically forget it ;-) – Przemyslaw Szufel Aug 24 '20 at 08:56
0

I guess you don't mean it like this, but one can also interpret the question in the way that y is a sparse vector in the Julia sense, and you want to construct a sparse diagonal matrix out of it. In that case you can do the following:

julia> y = sprand(10, 0.2)
10-element SparseVector{Float64,Int64} with 2 stored entries:
  [4 ]  =  0.389682
  [5 ]  =  0.232429

julia> I, V = findnz(y)
([4, 5], [0.3896822408908356, 0.2324294021548845])

julia> sparse(I, I, V)
5×5 SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [4, 4]  =  0.389682
  [5, 5]  =  0.232429

Unfortunately, spdiagm does not preserve structural zeros for a sparse input:

julia> spdiagm(0 => y)
10×10 SparseMatrixCSC{Float64,Int64} with 10 stored entries:
  [1 ,  1]  =  0.0
  [2 ,  2]  =  0.0
  [3 ,  3]  =  0.0
  [4 ,  4]  =  0.389682
  [5 ,  5]  =  0.232429
  [6 ,  6]  =  0.0
  [7 ,  7]  =  0.0
  [8 ,  8]  =  0.0
  [9 ,  9]  =  0.0
  [10, 10]  =  0.0

I don't know whether this is intentional, but I filed an issue about this behaviour.

phipsgabler
  • 20,535
  • 4
  • 40
  • 60
  • Note that `Diagonal` or `BandedMatrix` can have more efficient memory storage (because they explicitly know the data structure and support only this data structure). Hence while `SparseArrays` is far more flexible this also must come with significant performance cost. – Przemyslaw Szufel Aug 24 '20 at 08:58