# PSD sparse matrix generation: (almost) negative eigenvalues

**URL:** <https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851>\
**Category:** Optimization (Mathematical)\
**Created:** [June 26, 2023, 5:10pm UTC](https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851 "2023-06-26T17:10:33Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![selvi-aras](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/selvi-aras/32/50637_2.png) [@selvi-aras](https://discourse.julialang.org/u/selvi-aras)\
**Post date:** [June 26, 2023, 5:10pm UTC](https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851/1 "2023-06-26T17:10:33Z")

</div>

This question is not fully related to optimization (though I am using it for optimization); but a quick search on the Julia forum hinted me that such questions go to the optimization tag.

I can generate PSD matrices with:

```julia
n = 10
Q = randn(n,n); # random iid sample Gaussian elements
Q = Q'*Q; #PSD matrix generation
Q = (Q + Q')/2; # to ensure symmetry (sometimes numerical issues break it)

```

Now, I add a sparsification:

```julia
n = 10
Q = randn(n,n); 
Q = Q.*(rand(n,n) .> 0.5) # ******new******
Q = Q'*Q; 
Q = (Q + Q')/2; 

```

Sometimes the minimum eigenvalue of such generated `Q` is something like -10^{-8} and my solvers think this is not PSD. Typically adding a small diagonal matrix of 10^{-6} would be sufficient, but in this example, that does not save some instances. I wanted to check if there is any systematic sparse PSD matrix generation.

**Edit** : I fixed the issue by controlling the eigenvalues:

```julia
e_vals = 9 .* rand(n) .+ 1; #between 1-10.
A, _ = qr(randn(n, n) .* (rand(n,n) .> 0.75)); 
D = Matrix(Diagonal(e_vals)); Q = round.(A*D*A', digits = 8)

```

but still wondering if my initial example would also work.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [June 27, 2023, 3:09pm UTC](https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851/2 "2023-06-27T15:09:40Z")

</div>

The fixed version is generating well conditioned positive definite matrices, not rank deficient semidefinite matrices. That seems like it wouldn’t give you what you want if you really want semidefinite matrices that could be rank deficient.

In the second example, the sparsification is introducing rank deficiency in Q in some cases, depending on what gets set to zero. So you do get semidefinite matrices sometimes. But what it means for a matrix to be semidefinite in the presence of rounding is a bit tricky.

Even if your matrix really is truly rank deficient semidefinite, which is hard to do in this way with rounding unless the sparsity structure of Q^TQ enforces rank deficiency, an eigenvalue computation will at best give you the eigenvalues of Q+E\_1 where \|E\_1\| can be bounded in terms of the unit roundoff and \|Q\|. So you can have a PSD matrix Q and still compute a small negative eigenvalue of Q+E\_1. (Although -10^{-8} seems too large in magnitude unless you did some 8 digit rounding as in your fixed example; -10^{-16} is closer to what I might expect in double precision). And even if your matrix is PSD and you compute nonnegative eigenvalues, if it’s nearly rank deficient, whatever code you are using afterward might also solve some problem for a different perturbed matrix Q+E\_2 that is not positive semidefinite. So you will still run into trouble potentially. In practice, when rounding is involved, being to close to an indefinite matrix can be just as bad as being indefinite.

Of course if you are fine with generating well conditioned positive definite matrices, your fixed example works. Otherwise, modifying Q with a diagonal shift can work in some contexts. Or maybe the code you are feeding into needs to be more robust in handling matrices with small negative eigenvalues. There’s no generally applicable rule for this sort of problem. It depends on the details of what you are trying to do with Q.

---

<div class="post-metadata">

**Author:** ![fph](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fph/32/17159_2.png) [@fph](https://discourse.julialang.org/u/fph)\
**Post date:** [June 27, 2023, 3:49pm UTC](https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851/3 "2023-06-27T15:49:29Z")

</div>

> [@selvi-aras](#):
>
> ```julia
> Q = Q'*Q; #PSD matrix generation
> Q = (Q + Q')/2; # to ensure symmetry (sometimes numerical issues break it)
> 
> ```

Wait, what? I’m not a Blas/Lapack pro, but to me this seems like a bug in Julia base. The operation `Q' * Q` should map to a BLAS `syrk`, which is guaranteed to preserve symmetry because it computes off-diagonal values only once and then mirrors them. It works like that in Matlab.

If the results do not come up symmetric, this suggests that it’s mapping to `gemm` or `symm` instead, and it’s doing twice as many flops as necessary.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [June 27, 2023, 4:31pm UTC](https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851/4 "2023-06-27T16:31:03Z")

</div>

The debugger shows that it eventually calls `syrk`. That’s nice, but I’m not sure I would have assumed that it would do that or something comparable now and forever unless I saw it documented somewhere. So I might still be tempted to do something like that if I was really depending on exact symmetry.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 27, 2023, 5:24pm UTC](https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851/5 "2023-06-27T17:24:47Z")

</div>

> [@selvi-aras](#):
>
> ```julia
> Q = randn(n,n); 
> Q = Q.*(rand(n,n) .> 0.5) # ******new******
> 
> ```

As an aside, note that this is statistically equivalent to

```julia
Q = sprand(n, n, 0.5, randn)

```

which is much more efficient _if_ you go to larger matrices and lower probabilities (so that sparse-matrix data structures become beneficial). Of course, in this case it’s even more likely to have zero eigenvalues (which roundoff errors may make slightly negative).

---

<div class="post-metadata">

**Author:** ![selvi-aras](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/selvi-aras/32/50637_2.png) [@selvi-aras](https://discourse.julialang.org/u/selvi-aras)\
**Post date:** [June 28, 2023, 12:13am UTC](https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851/6 "2023-06-28T00:13:43Z")

</div>

I remember very specifically yesterday that without the sum / 2 trick the symmetry was not recognized. I will have a look again – I hope I am not misleading anyone.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [June 28, 2023, 2:52am UTC](https://discourse.julialang.org/t/psd-sparse-matrix-generation-almost-negative-eigenvalues/100851/7 "2023-06-28T02:52:43Z")

</div>

> [@fph](#):
>
> The operation `Q' * Q` should map to a BLAS `syrk`, which is guaranteed to preserve symmetry because it computes off-diagonal values only once and then mirrors them.

No, `syrk` does not do such a mirror/copy (at least in the reference implementation), although in this context (`Q'*Q`) a Julia wrapper does. The BLAS convention is that only either upper or lower triangular elements of the result are defined, so no copy is needed. Whereas symmetry is unsafely defined by contract in C or Fortran BLAS, it could be enforced in Julia by the type system, viz. returning `Symmetric` in this context, also avoiding the copy. @stevengj do you know if that option is unwise or has been discussed somewhere?