# MATLABs \`sprandsym\` in Julia

**URL:** <https://discourse.julialang.org/t/matlabs-sprandsym-in-julia/72251>\
**Category:** General Usage\
**Tags:** question, matlab, sparse\
**Created:** [November 29, 2021, 6:47pm UTC](https://discourse.julialang.org/t/matlabs-sprandsym-in-julia/72251 "2021-11-29T18:47:52Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [November 29, 2021, 6:47pm UTC](https://discourse.julialang.org/t/matlabs-sprandsym-in-julia/72251/1 "2021-11-29T18:47:52Z")

</div>

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)`](https://www.mathworks.com/help/matlab/ref/sprandsym.html), 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.

```julia
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.)

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [November 29, 2021, 6:50pm UTC](https://discourse.julialang.org/t/matlabs-sprandsym-in-julia/72251/2 "2021-11-29T18:50:39Z")

</div>

FWIW, another (orthogonal) idea based on [How to generate random symmetric positive definite matrices using MATLAB? - Mathematics Stack Exchange](https://math.stackexchange.com/a/358092) that I tried was

```julia
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.