Spdiag Matlab function in Julia

General Usage
1

Is it possible to create an upper triangular matrix in Julia as fast as in Matlab, given the vectors that occupy “diagonal” positions? To be precise, spdiags allows to locate a vector on the diagonal (index=0) or the first diagonal on the right of the matrix diagonal (index=1), this way creating an upper diagonal matrix, and so on and so forth. The only information needed are the vector itself and the position (0,1,2,3,…N-1, -1,-2,-3,…) where negative numbers are positions in a lower triangular matrix. No need to provide the entire vector of indeces!

2

Yes, you are looking for spdiagm in Julia. You need to add using SparseArrays. For more information see: Sparse Arrays · The Julia Language . Unfortunately, it creates a square matrix, so you need to cut off the lines or columns that you don’t want. Here one example:

MVtmp = spdiagm(2 => ones(4), 0 => -ones(4));
Matrix(MVtmp)

Output:
6×6 Array{Float64,2}:
 -1.0   0.0   1.0   0.0  0.0  0.0
  0.0  -1.0   0.0   1.0  0.0  0.0
  0.0   0.0  -1.0   0.0  1.0  0.0
  0.0   0.0   0.0  -1.0  0.0  1.0
  0.0   0.0   0.0   0.0  0.0  0.0
  0.0   0.0   0.0   0.0  0.0  0.0

but if you want a 4x6 matrix, you need to do:

MV = MVtmp[1:4,:];
Matrix(MV)

Output:
4×6 Array{Float64,2}:
 -1.0   0.0   1.0   0.0  0.0  0.0
  0.0  -1.0   0.0   1.0  0.0  0.0
  0.0   0.0  -1.0   0.0  1.0  0.0
  0.0   0.0   0.0  -1.0  0.0  1.0

I’ve used the function Matrix() to visualise the sparse matrix.

3 Likes
3

This was fixed in Julia 1.3 — spdiagm now allows you to pass optional size arguments to create non-square matrices:

julia> Matrix(spdiagm(4, 6, 2 => ones(4), 0 => -ones(4)))
4×6 Array{Float64,2}:
 -1.0   0.0   1.0   0.0  0.0  0.0
  0.0  -1.0   0.0   1.0  0.0  0.0
  0.0   0.0  -1.0   0.0  1.0  0.0
  0.0   0.0   0.0  -1.0  0.0  1.0
6 Likes
4

Wonderful!

5

What if I want to declare a matrix that’s 1000x1000 with a hundred populated diagonals? spdiags in Matlab can take a dense array that stores the diagonal entries, a vector of offsets, and matrix dimension paramters to build n-diagonal matrix using: A = spdiags(Bin,d,m,n);

Is there a comparable syntax for this?

6

You would call spdiagm from the SparseArrays stdlib, e.g. you can emulate the Matlab syntax A = spdiags(Bin,d,1000,1000) with:

spdiagm(1000,1000, [d[i] => Bin[1:1000-abs(d[i]), i] for i = 1:length(d)]...)

Note that Matlab’s Bin is a bit odd because Matlab only uses part of each column Bin[:,i] depending on the diagonal d[i], whereas Julia expects a vector of the correct length for the corresponding diagonal.

The code is a bit simpler if you have an array of arrays B, where each element B[i] is a vector (of the correct length) of the entries on diagonal i:

spdiagm(1000,1000, (d .=> B)...)

Here, I’m using the broadcasting “dot” operator .=> version of the pair constructor =>.

For example:

julia> using SparseArrays

julia> d = [7,3,0,-2,-8];

julia> B = [rand(1000 - abs(d)) for d in d];

julia> A = spdiagm(1000,1000, (d .=> B)...)
1000×1000 SparseMatrixCSC{Float64, Int64} with 4980 stored entries:
⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦

(The last output is the new sparse-matrix display format in Julia 1.6.)

You should also look into the BandedMatrices.jl package, which is likely to be much more efficient than generic sparse matrices in this case, especially since your matrix is not particularly sparse (10% sparsity).

3 Likes
7

I love the .=> notation, though splatting still confuses me. I was going to kludge a packed array to I, J, V converter, but this is much better. Thanks!

Just for those wanting to translate Matlab code directly, the syntax for your first solution should be:
spdiagm(1000,1000, [d[i] => Bin[1+maximum([0,d[i]]):1000+minimum([0,d[i]]), i] for i = 1:length(d)]...)
This takes care of the way Matlab wants the zero padding of the packed array to work.

1 Like
8

Better to just do max(0,d[i]) and min(0,d[i]).

(This is a Matlab-ism and Julia anti-pattern — allocating an array of 2 elements just to take the maximum. Matlab (and Numpy and R) really train programmers in a distorted style, unfortunately.)

2 Likes