# MATLABs `sprandsym` in Julia

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.

(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.)
``````A = sprandn(n, n, p)