MATLABs `sprandsym` in Julia

General Usage
, ,
1

Hello everyone!

I need a function that generates a random sparse matrix of size (n,n) with the following properties:

  • issymmetric
  • isposdef
  • has approximately p * n^2 non-zero elements (for given p)
  • has a condition number of exactly c (for given c)

In MATLAB, there is sprandsym(n,density,rc,kind), which gives me such a matrix. However, I couldn’t find a Julia pendant, which is why I started to write one.

I focused on the idea mentioned in the docs of MATLABs sprandsym (kind=1), that is to start with a positive definite diagonal matrix with the desired condition number and then apply Jacobi rotations to generate non-zero off-diagonal elements. However, I fail to see how to efficiently get the desired density. Currently, my best attempt is the following.

function sprandsymposdef(n, p, c=1)
    # check inputs
    0 ≤ p ≤ 1 || throw(ArgumentError("density must be ≥ 0 and ≤ 1"))
    c ≥ 1 || throw(ArgumentError("condition number must be ≥ 1"))
    n ≥ 1 || throw(ArgumentError("n must be ≥ 1"))

    # enforce condition number
    A = spdiagm(n, n, range(1, c; length=n))

    # apply random Jacobi rotations
    #     - creates non-zero off diagonal elements
    #     - preserves eigvals, thus also posdef property and cond
    # TODO: speed up / avoid allocs by updating the elements of A
    #       directly (i.e. don't construct Jacobi matrix explictly.)
    while density(A) ≤ p
        a, b = (rand(1:n), rand(1:n))
        J = jacobi_rot_mat(n, a, b)
        A = J' * A * J
    end

    # enforce symmetry, i.e. iron out numerical inaccuracies
    A = (A + A') .* 0.5
    return A
end

function jacobi_rot_mat(n, p, q)
    c = rand()
    s = sqrt(1 - c^2)
    J = sparse(I(n) * 1.0)
    J[p, p] = c
    J[q, q] = c
    J[q, p] = s
    J[p, q] = -s
    return J
end

While this works, often tens or hundreds of thousands of Jacobi rotations are necessary and, even worse, this number strongly depends on the values of p and n such that the performance varies a lot for different inputs. (A reasonable input could be sprandsymposdef(10000, 0.001, 2) or sprandsymposdef(10000, 0.01, 2) to really see the issue when trying to reach a higher density.)

So, I really hope that you can help me improve my strategy or perhaps even figure out an entirely different approach. Any kind of help is very much appreciated.

Thanks in advance!

(Note that, obviously, I can speed up the code above by avoiding the construction of the Jacobi rotation matrix and instead updating A in-place etc. However, I’m really hoping for an algorithmic improvement / suggestion before I start optimizing the computation itself.)

2

FWIW, another (orthogonal) idea based on How to generate random symmetric positive definite matrices using MATLAB? - Mathematics Stack Exchange that I tried was

A = sprandn(n, n, p)
A = sparse(Symmetric(A))
A = A + n * I(n)

which checks many of the boxes (including the approximate density). However, here I don’t know how to reliably control the condition number.

2 Likes