# Issue: PDMats.jl matrices losing positive-definiteness with Optim.jl autodiff

**URL:** https://discourse.julialang.org/t/issue-pdmats-jl-matrices-losing-positive-definiteness-with-optim-jl-autodiff/133421
**Category:** General Usage
**Tags:** question
**Created:** [October 24, 2025, 10:50pm UTC](https://discourse.julialang.org/t/issue-pdmats-jl-matrices-losing-positive-definiteness-with-optim-jl-autodiff/133421 "2025-10-24T22:50:52Z")
**Posts on this page:** 1
**Showing post:** 4

<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: [October 25, 2025, 11:44am UTC](https://discourse.julialang.org/t/issue-pdmats-jl-matrices-losing-positive-definiteness-with-optim-jl-autodiff/133421/4 "2025-10-25T11:44:19Z")

</div>

> [@kellertuer](#):
>
> An advantage of this approach is, that you suddenly can use **unconstraint** optimization like a L-BFGS on that manifold. The replacement for + (and -) is the exponential (and logarithmic) map, they are a bit more expensive than +/-, so there is no free lunch.

You could also just use a change of variables to do nonlinear optimization over semidefinite matrices (which you have to include in “SPD” for optimization since the set of strictly positive-definite matrices is not closed). That is, instead of optimizing f(A) over “SPD” matrices A, you optimize f(X^T X) over _unconstrained_ m \times m matrices X, or better yet `f(Hermitian(X'X))` to ensure that it is exactly Hermitian despite roundoff errors.

This has no loss of generality, allows you to use _any_ nonlinear optimization algorithm, and is easy. (The downside is that you might lose some special structure in f that would otherwise let you use something specialized like an SDP solver, but I’m guessing you don’t have a straightforward SDP or you would be using such a solver already.)

You can also, of course, use specialized methods for semidefinite constraints. And, of course, it depends a lot on what you want to _do_ with the SPD matrix in your function. If you want to solve linear systems, for example, it would be more efficient to optimize over the vector space of lower- or upper-triangular matrices X, since then you have the Cholesky factorization for free.

Note, however, that X^T X may be slightly indefinite due to roundoff errors, so your f(X^T X) computation needs to be robust to such things. This is a general issue to be aware of when working with SDP matrices; how to ensure robustness depends on what precise calculation you’re performing, but a general idea is to treat slightly negative eigenvalues as zeros. Alternatively you can add some regularization like replacing f(A) with f(A + \varepsilon \Vert A \Vert\_F I / m) for a small \varepsilon proportional to the precision. This has come up many times on this forum, e.g.:

- [Numerical inaccuracies when constructing correlation matrices - #6 by stevengj](https://discourse.julialang.org/t/numerical-inaccuracies-when-constructing-correlation-matrices/85177/6)
- [Random draws of multivariate normal with positive semi-definite covariance matrix - #2 by stevengj](https://discourse.julialang.org/t/random-draws-of-multivariate-normal-with-positive-semi-definite-covariance-matrix/78493/2)
- [Cholesky Decomposition of a Sparse Symmetric Positive Semidefinite (SPSD) Singular Matrix - #8 by stevengj](https://discourse.julialang.org/t/cholesky-decomposition-of-a-sparse-symmetric-positive-semidefinite-spsd-singular-matrix/119682/8)
- [Cholesky decomposition of low-rank positive-semidefinite matrix](https://discourse.julialang.org/t/cholesky-decomposition-of-low-rank-positive-semidefinite-matrix/70397)
- [What is the fastest way to check if a hermitian matrix is positive semi-definite? - #8 by mstewart](https://discourse.julialang.org/t/what-is-the-fastest-way-to-check-if-a-hermitian-matrix-is-positive-semi-definite/110769/8)

(But this may not be an issue if your optimization problem does not favor matrices that are close to singular, i.e. if the SDP constraint is inactive at the optimum — in that case, as long as you start with a random X so that X^T X is far from singular, probably the gradient steps will keep things that way.)

---

_[View the full topic](https://discourse.julialang.org/t/issue-pdmats-jl-matrices-losing-positive-definiteness-with-optim-jl-autodiff/133421)._