# How to restrict a matrix to be positive definite

**URL:** <https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318>\
**Category:** General Usage\
**Tags:** question, math\
**Created:** [February 26, 2025, 12:50am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318 "2025-02-26T00:50:17Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [February 26, 2025, 12:50am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/1 "2025-02-26T00:50:17Z")

</div>

I am estimating a model. In the process of estimation I have to convert a fixed independent and normally distributed vector to a correlated multivariate normal vector. The way I am doing that is using this:

\bf{x} \sim N(\bf{0}, I) \Rightarrow S \bf{x} \sim N(0,SS').

That means that if I can obtain a matrix \bf{S} such that \bf{S}\bf{S}' is positive semi definite, I can premultiply vector \bf{x} by \bf{S} to obtain a correlated normal vector. The problem I’m facing is that I am parameterizing matrix \Sigma = \bf{S} \bf{S}'. Say it is 3x3:

```julia
Σ = [s11 s12 s13;
     s12 s22 s23;
     s13 s23 s33]

```

Those are the parameters I am supplying the optimizer. I am obtaining \bf{S} using a Cholesky decomposition.

How can I impose restrictions on the components of `Σ` such that it is positive semi definite (and therefore the Cholesky decomposition exists)?

I thought about parameterizing \bf{S} instead, and then constructing `Σ = S * transpose(S)`. The problem in that case is that `Σ` ends up having all positive values and I want to allow for some negative components in it (allowing for negative covariance between components of \bf{x}).

Here is a simple example:

```julia
using LinearAlgebra

x = randn(3, 100)

function obj(s; x)
    s11 = s[1]
    s12 = s[2]
    s13 = s[3]
    s22 = s[4]
    s23 = s[5]
    s33 = s[6]

    Σ = [s11 s12 s13;
         s12 s22 s23;
         s13 s23 s33]
    S = cholesky(Σ).L
    y = S * x 
    return f(y)
end

f(y) = sum(y)

```

Again, I want to be able to supply function `obj` to an optimizer, which means that it should be able to accept any vector `s`. What are restrictions I can place on `Σ` to ensure it is positive semi definite? For example, the main diagonal elements must be non-negative, which can be achieved with an `exp` transform.

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [February 26, 2025, 12:55am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/2 "2025-02-26T00:55:46Z")

</div>

Use the function `isposdef` after `using LinearAlgebra`?

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [February 26, 2025, 12:57am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/3 "2025-02-26T00:57:27Z")

</div>

That only checks for psd, it does not restrict `Σ` to be one.

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [February 26, 2025, 1:03am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/4 "2025-02-26T01:03:06Z")

</div>

What do you mean by restricting it to be one then?

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [February 26, 2025, 1:33am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/5 "2025-02-26T01:33:44Z")

</div>

What I want is to make function `obj` able to accept any inputs, transform them, and make a psd matrix from them. That way, I can safely supply function `obj` to an optimizer without having to worry whether the function will error because it can’t construct psd `Σ`.

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [February 26, 2025, 1:37am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/6 "2025-02-26T01:37:31Z")

</div>

In that case, something like starting with the Σ matrix that you have. Then,

Σ = Symmetric( Σ+Σ’) / 2

to make it symmetric. Then compute its minimum eigenvalue λ.

Then add (c-λ) \* I to Σ

where c is the desired smallest eigenvalue (some positive number)

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [February 26, 2025, 1:47am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/7 "2025-02-26T01:47:25Z")

</div>

Thank you. Does that method “span” all positive psd matrices? That is, can I construct every possible psd matrix by varying c and Σ?

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [February 26, 2025, 2:01am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/8 "2025-02-26T02:01:56Z")

</div>

Yes

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [February 26, 2025, 3:29am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/9 "2025-02-26T03:29:33Z")

</div>

Since c is a parameter of my choice, do I also include it with the rest, so that the optimizer can control it?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 26, 2025, 4:20am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/10 "2025-02-26T04:20:18Z")

</div>

> I thought about parameterizing instead, and then constructing `Σ = S * transpose(S)`. The problem in that case is that `Σ` ends up having all positive values and I want to allow for some negative components in it (allowing for negative covariance between components of ).

I think this is probably the better approach. That said, I’m not sure why you would be seeing `Σ` ending up having all positive values. Were you making `S` be symmetric rather than triangular?

```julia
julia> S = [1 -1
            0 1];

julia> S*S'
2×2 Matrix{Int64}:
  2 -1
 -1 1

```

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [February 26, 2025, 5:51am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/11 "2025-02-26T05:51:54Z")

</div>

Ah yes, I was making that mistake!

---

<div class="post-metadata">

**Author:** ![Alexander\_Knudson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alexander_knudson/32/215656_2.png) [@Alexander\_Knudson](https://discourse.julialang.org/u/Alexander_Knudson)\
**Post date:** [February 26, 2025, 3:07pm UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/12 "2025-02-26T15:07:09Z")

</div>

I’m not exactly sure how to solve your problem, but there’s the NearestCorrelationMatrix.jl library that can find a positive definite correlation matrix from an input matrix. Can you compute the the covariance matrix from the correlation matrix?

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [February 26, 2025, 4:52pm UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/13 "2025-02-26T16:52:35Z")

</div>

For optimization over constrained spaces, I usually reach for [TransformVariables.jl](https://github.com/tpapp/TransformVariables.jl) providing transformations onto unconstrained vectors. In particular, for covariance matrices check out [corr\_cholesky\_factor](https://www.tamaspapp.eu/TransformVariables.jl/stable/#TransformVariables.CorrCholeskyFactor).  
In case you want to get really fancy, you can also use [Manopt.jl](https://manoptjl.org/stable/) and optimize over the manifold of [symmetric positive definite matrices](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/symmetricpositivedefinite/) directly.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [March 2, 2025, 6:03am UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/14 "2025-03-02T06:03:31Z")

</div>

I was checking `TransformVariables.jl`, the author says that the densities need to be adjusted. Why is that necessary? What I had in mind to do is to employ the delta method to obtain the correct variances of the transformed parameters (those in \Sigma) from the obtained estimates (those in S).

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [March 5, 2025, 5:24pm UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/15 "2025-03-05T17:24:20Z")

</div>

When you change your parametrization, e.g., mapping from an unconstrained vector x \in \mathbb{R}^n to some constrained object \theta \in \Theta, the density changes accordingly. Let g: x \to \theta, then

p(x) = p(\theta) \left| \frac{d\theta}{dx} \right|

where \left| \frac{d\theta}{dx} \right| denotes the log Jacobian of g.  
E.g., in Bayesian estimation this is important when the prior is defined on \theta, but the sampling is done on x. When optimizing the likelihood this is usually not required as the maximum is not effected by the coordinate transform, e.g., \mathrm{argmax}\_{\theta} p(D | \theta) = g\left( \mathrm{argmax}\_{x} p(D | g(x))\right).  
Don’t know how/if that applies to the delta method … hope it helps anyways.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [March 5, 2025, 10:35pm UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/16 "2025-03-05T22:35:46Z")

</div>

Thank you! Is there any name for that? Can you provide some references where I can read more about it?

What I had in mind to do was maximize the likelihood with respect to unconstrained vector, obtain the variance covariance of those parameters, and then use the delta method to obtain the variance covariance for the transformed parameters:

\hat{x} = \arg \max\_x \ell(p(x)).

The variance covariance of \hat{x} is

var(\hat{x}) = -\left(\frac{\partial^2{\ell}}{\partial x^2}(\hat{x})\right)^{-1}.

Ultimately, I am interested in the variance covariance of the transformed parameters \hat{\theta} = p(\hat{x}). For that, I apply the delta method:

\sqrt{n} (p(\hat{x}) - p(x)) \sim N(0, h'var(\hat{x})h),

where h = \frac{\partial}{\partial x}p(x)'.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [March 6, 2025, 6:00pm UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/17 "2025-03-06T18:00:31Z")

</div>

In the context of Bayesian modeling, the [Stan manual](https://mc-stan.org/docs/stan-users-guide/reparameterization.html) discusses reparametrization and change of variables.  
If I get your case correctly, you maximize the likelihood on the unconstrained space, approximate the uncertainty there via the Laplace approximation, i.e., using the Hessian as inverse variance, and then use another approximation – the delta method – in order to obtain the variance on the constrained space? Here, I would probably skip the delta method altogether and simply use the Laplace approximation on the constraint space directly, i.e., compute \frac{\partial^2 \mathcal{l}}{\partial \theta}(\hat{\theta}) directly. In any case, you won’t need to correct for the change of variable as the transformation is applied to the point estimates \hat{x} \to \hat{\theta} only.

---

<div class="post-metadata">

**Author:** ![Philippe\_Maincon1](https://avatars.discourse-cdn.com/v4/letter/p/ec9cab/32.png) [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)\
**Post date:** [March 6, 2025, 6:38pm UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/18 "2025-03-06T18:38:34Z")

</div>

Any symmetric positive definite (SPD) matrix has a unique Cholesky decomposition  
`M=L*L'`  
where `L` is lower triangular. Hence the space of `n` by `n` SPD matrices has dimension `n*(n+1)/2`, and a way to constrain an algorithm to only explore all such matrices and only them is to use the upper triangular coefficients of `L` as unknowns.

---

<div class="post-metadata">

**Author:** ![npn0010](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/npn0010/32/44420_2.png) [@npn0010](https://discourse.julialang.org/u/npn0010)\
**Post date:** [September 5, 2025, 4:38pm UTC](https://discourse.julialang.org/t/how-to-restrict-a-matrix-to-be-positive-definite/126318/19 "2025-09-05T16:38:59Z")

</div>

Have you considered using eigendecomposition to parameterize \Sigma? I.e., you can let \Sigma=Q\Lambda Q^\top where \Lambda is a diagonal matrix of eigenvalues (set to be positive definite to ensure \Sigma is positive definite) and where Q is an orthogonal matrix parameterized in some particular way. See [this arxiv paper](https://arxiv.org/pdf/2504.16328) for an example of this.