# Exploring Quadratic forms with Julia

**URL:** <https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138>\
**Category:** Teaching & Outreach\
**Created:** [September 22, 2023, 9:53am UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138 "2023-09-22T09:53:18Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![VinodV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinodv/32/32733_2.png) [@VinodV](https://discourse.julialang.org/u/VinodV)\
**Post date:** [September 22, 2023, 9:53am UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/1 "2023-09-22T09:53:18Z")

</div>

Hi friends,  
I have published a blog post on [Exploring Quadratic forms with Julia](https://forem.julialang.org/vinodv/exploring-quadratic-forms-with-julia-k6e)  
Just tried to make my understanding on Quadratic forms clearer using Julia.  
Please give your feedback on improvements/efficient code etc.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [September 23, 2023, 11:05pm UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/2 "2023-09-23T23:05:59Z")

</div>

A few comments:  
For convenience, I copy&paste the code from Forem here:

```plaintext
using LinearAlgebra
function getDefiniteness(A)

    @assert isequal(A,A') "Matrix is not symmetric"

    if true == isposdef(A)
        return "Matrix is positive definite"
    end

    x,y = eigmin(A),eigmax(A)

    if iszero(x)
        return "Matrix is positive semi definite"
    end

    if y > 0.0 && x < 0.0
        return "Matrix is indefinite"
    end

    if iszero(y)
        return "Matrix is negative semi definite"
    end

    if all([x, y] .< 0)
        return "Matrix is negative definite"
    end

end

A = rand(5,5)
A = (A+A')/2

A = [-1 0;0 -1.0]
getDefiniteness(A)

```

First, you do not have to write

```julia
if true == isposdef(A)

```

Instead, you can write directly

```julia
if isposdef(A)

```

Second, since you compute eigenvalues for your matrix in order to distinguish between definiteness and indefiniteness, you can skip the the `isposdef()` computation altogether to save some computational effort (should that be an issue in a particular situation). Once you have the eigenvalues, you can use them for checking indefiniteness too.

Third, in finite-precision arithmetics, it rarely makes any sense to check if something is exactly zero as you do with

```plaintext
if iszero(x)

```

Consider the following example

```julia
A = [1e-100 0;0 1.0]

```

The matrix is classified as positive definite by your function

```plaintext
julia> getDefiniteness(A)
"Matrix is positive definite"

```

while it is effectively/practically singular (positive semidefinite) for all intents and purposes.

Finally, although I am not familiar with the implementation details of `eigmin`, by measuring the execution times and seeing that they are about the same as those for `eigvals`, I suspect that all eigenvalues are computed even if you only request the smallest one. I conclude that instead of calling `eigmin(A)` and then `eigmax()`, it is perhaps more efficient to call `eigvals` just once…

---

<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:** [September 24, 2023, 1:46am UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/3 "2023-09-24T01:46:28Z")

</div>

> [@zdenek\_hurak](#):
>
> ` @assert isequal(A,A') "Matrix is not symmetric"`

You can just use `ishermitian(A)`. And after that I would set `H = Hermitian(A)` and use `H` instead of `A` (so that Julia will know to use more efficient algorithms).

And normally you would not use `@assert` for this — you should assume that `@assert` is something that might be removed in production code, and it should only be used to check for bugs. Throw an exception instead, e.g.

```julia
ishermitian(A) || throw(ArgumentError("matrix must be Hermitian"))

```

> [@zdenek\_hurak](#):
>
> `if true == isposdef(A)`

It is more idiomatic to simply do `if isposdef(A)` (or `isposdef(H)` if you use the Hermitian wrapper).

> [@zdenek\_hurak](#):
>
> `if all([x, y] .< 0)`

This allocates two arrays (`[x, y]` is one array, and `[x, y] .< 0` is another array), just to check if two numbers are negative. Putting everything into an array, even if it is just two scalars, is kind of a Matlab anti-pattern. Just check `if x < 0 && y < 0`.

> [@zdenek\_hurak](#):
>
> Second, since you compute eigenvalues for your matrix in order to distinguish between definiteness and indefiniteness, you can skip the the `isposdef()` computation altogether to save some computational effort

I would do the opposite: if you want to check whether a matrix is negative definite, it’s probably more efficient to call `isposdef(-A)` or `isposdef(-H)` than to look at eigenvalues. The reason is that you can check positive-definiteness using only a Cholesky factorization, which is _much_ cheaper than computing eigenvalues.

For example, with a 1000 \times 1000 negative-definite matrix:

```julia
julia> A = randn(1000,1000); A = -A'A; H = Hermitian(A);

julia> @btime isposdef($H) || isposdef(-$H);
  5.998 ms (6 allocations: 22.89 MiB)

julia> @btime eigvals($H);
  42.406 ms (11 allocations: 7.99 MiB)

julia> @btime (eigmin($H), eigmax($H));
  67.561 ms (22 allocations: 15.96 MiB)

julia> @btime eigmin($H);
  31.952 ms (11 allocations: 7.98 MiB)

```

Notice that checking _both_ `isposdef(H)` and `isposdef(-H)` is 7x faster than computing `eigvals(H)`.

> [@zdenek\_hurak](#):
>
> Finally, although I am not familiar with the implementation details of `eigmin`, by measuring the execution times and seeing that they are about the same as those for `eigvals`, I suspect that all eigenvalues are computed even if you only request the smallest one.

No (at least not for a Hermitian matrix). Notice that for my example above, `eigmin(H)` is about 25% faster than `eigvals(H)`. However, it is still slower to call _both_ `eigmin` and `eigmax`, since they don’t share any computations.

In principle, it should be possible to compute _both_ `eigmin` and `eigmax` more efficiently than computing all the eigenvalues, but you need to share the Hessenberg/tridiagonal reduction calculation (since that is the first step of any eigenvalue calculation, and is a big part of the cost):

```julia
function eigminmax(H::Hermitian)
    T = hessenberg(H).H # real SymTridiagonal reduction
    return eigmin(T), eigmax(T)
end

```

which gives

```julia
julia> @btime eigminmax($H);
  32.574 ms (65 allocations: 8.27 MiB)

```

which is almost twice as fast as calling `(eigmin(H), eigmax(H))`, about the same speed as `eigmin(H)`, and 25% faster than `eigvals(H)`. However, it is still much slower than `isposdef(H)` (which computes no eigenvalues at all).

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [September 24, 2023, 9:24am UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/4 "2023-09-24T09:24:59Z")

</div>

I perfectly agree that Cholesky-decomposition-based `isposdef(H)` is the way to go for checking positive definiteness of `H` (and negative definiteness of `-H`). But in the OP’s code from Forem that I pasted here, it is not only positive definiteness and negative definiteness, but also positive semidefiniteness, negative semidefiniteness and indefiniteness that must be determined. I think this imposes the need to compute eigenvalues, doesn’t it? I do not know how I can distinguish between (strict) positive semidefiniteness and nondefiniteness when detecting the failure to compute the Cholesky decomposition. Then my recommendation was simply: once you are forced to compute the eigenvalues, you can use them to check for pd-ness and nd-ness too.

---

<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:** [September 25, 2023, 8:36am UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/5 "2023-09-25T08:36:28Z")

</div>

> [@zdenek\_hurak](#):
>
> positive definiteness and negative definiteness, but also positive semidefiniteness, negative semidefiniteness and indefiniteness

You don’t need eigenvalues. If it’s not positive-definite or negative-definite, you should assume it’s indefinite.

You can’t check for semi-definiteness numerically, because roundoff errors usually make it impossible to distinguish a zero eigenvalue from one that is slightly positive or slightly negative.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [September 25, 2023, 9:13am UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/6 "2023-09-25T09:13:17Z")

</div>

> [@stevengj](#):
>
> > [@zdenek\_hurak](#):
> >
> > `@assert isequal(A,A') "Matrix is not symmetric"`
> 
> You can just use `ishermitian(A)`. And after that I would set `H = Hermitian(A)` and use `H` instead of `A` (so that Julia will know to use more efficient algorithms).

In addition, there are some pathological cases where `ishermitian` is **much** faster than `isequal`, e.g.

```julia
julia> using LinearAlgebra, SparseArrays, BenchmarkTools

julia> let A = sprandn(1000, 1000, 0.01);
           H = A + A'
           @btime ishermitian($H)
           @btime $H == $H'
       end;
  24.727 μs (1 allocation: 8.00 KiB)
  520.155 μs (0 allocations: 0 bytes)

```

I’ve been bit by this before with really large, very sparse matrices where I ended up spending the majority of my computation time doing what I thought was a very cheap verification that my array was Hermitian.

```julia
julia> let A = sprandn(100000, 100000, 0.001);
           H = A + A'
           @btime ishermitian($H)
           @btime $H == $H'
       end;
  97.176 ms (2 allocations: 781.36 KiB)
  1.974 s (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [September 25, 2023, 9:35am UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/7 "2023-09-25T09:35:55Z")

</div>

> [@stevengj](#):
>
> You can’t check for semi-definiteness numerically, because roundoff errors usually make it impossible to distinguish a zero eigenvalue from one that is slightly positive or slightly negative.

I understand the argument, but isn’t it the same situation as in determining the rank of a matrix? Say, if we want to determine if a matrix is singular or not. Strictly speaking, we would (effectively) never be able to detect a singular matrix numerically, wouldn’t we?

For example, is this matrix

```julia
julia> H = [-1e-20 0; 0 1]
2×2 Matrix{Float64}:
 -1.0e-20 0.0
  0.0 1.0

```

singular?

```julia
julia> rank(H)
1

```

For practical purposes, shall we regard it as indefinite or positive semi-definite?

---

<div class="post-metadata">

**Author:** ![VinodV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinodv/32/32733_2.png) [@VinodV](https://discourse.julialang.org/u/VinodV)\
**Post date:** [September 25, 2023, 2:53pm UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/8 "2023-09-25T14:53:44Z")

</div>

```julia
julia> @btime a,b = eigmin(A),eigmax(A)
  3.900 μs (23 allocations: 5.25 KiB)
(-1.4340659161787122, 5.636707875939123)

julia> @btime a = eigvals(A)
  1.880 μs (10 allocations: 2.59 KiB)
5-element Vector{Float64}:
 -1.4340659161787122
 -0.4391186490192085
  0.15889223870067884
  1.5551498036607247
  5.636707875939123

julia> @btime a = eigvals(A)
  1.880 μs (10 allocations: 2.59 KiB)
5-element Vector{Float64}:
 -1.4340659161787122
 -0.4391186490192085
  0.15889223870067884
  1.5551498036607247
  5.636707875939123

julia> @btime a,b = eigmin(A),eigmax(A)
  3.913 μs (23 allocations: 5.25 KiB)
(-1.4340659161787122, 5.636707875939123)

```

eigvals is efficient for a 5x5 matrix.

---

<div class="post-metadata">

**Author:** ![Syx\_Pek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/syx_pek/32/6364_2.png) [@Syx\_Pek](https://discourse.julialang.org/u/Syx_Pek)\
**Post date:** [September 25, 2023, 7:09pm UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/9 "2023-09-25T19:09:22Z")

</div>

`eigmin` and `eigmax` seems to call `eigvals` internally, so this makes sense (and should be true regardless of the size of matrix).

---

<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:** [September 25, 2023, 9:53pm UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/10 "2023-09-25T21:53:12Z")

</div>

> [@zdenek\_hurak](#):
>
> I understand the argument, but isn’t it the same situation as in determining the rank of a matrix? Say, if we want to determine if a matrix is singular or not. Strictly speaking, we would (effectively) never be able to detect a singular matrix numerically, wouldn’t we?

Right. You can only compute the “numerical rank” to some tolerance instead. See [this discussion](https://github.com/JuliaLang/julia/issues/48058#issuecomment-1371263157) and references therein.

In numerical computations, it’s more natural to talk about the **condition number** of a matrix, which indicates “how singular” a matrix is (and generalizes to non-square cases). Even if a matrix is non-singular, if it is _close enough_ to singular then you have to approach many calculations with care.

> [@zdenek\_hurak](#):
>
> `H = [-1e-20 0; 0 1]`
> 
> […] For practical purposes, shall we regard it as indefinite or positive semi-definite?

`cond(H) ≈ 1e20`, so for practical purposes you regard it as **ill-conditioned**.

---

<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:** [September 25, 2023, 10:02pm UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/11 "2023-09-25T22:02:04Z")

</div>

> [@Syx\_Pek](#):
>
> `eigmin` and `eigmax` seems to call `eigvals` internally, so this makes sense (and should be true regardless of the size of matrix).

They call `eigvals` but (for Hermitian/real-symmetric matrices) pass the [optional `irange` argument](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.eigvals) to compute only a subset of the eigenvalues. This saves some computation time, although the bulk of the time is spent on the Hessenberg/tridiagonal reduction as I [commented above](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/3) (so you ideally want to re-use that if you want _both_ `eigmin` and `eigmax`).

LAPACK isn’t very well optimized for tiny matrices like 5x5 so benchmarks are mostly measuring overheads; in principle you can do _much_ better for tiny matrices, e.g. using StaticArrays if the sizes are fixed.

---

<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:** [September 25, 2023, 11:20pm UTC](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138/12 "2023-09-25T23:20:57Z")

</div>

To some extent you can determine if a Hermitian matrix is close to semidefinite by computing the Cholesky decomposition with diagonal pivoting P^T H P = L L^T and using the largest magnitude pivot on the diagonal in each step. If your largest available pivot at step r is small, the matrix is correspondingly close to a semidefinite matrix of rank r. There is an analysis of the procedure in this [paper](https://eprints.maths.manchester.ac.uk/1193/1/high90c.pdf). The bounds on what you have to truncate aren’t as strong as with eigenvalues, but it works well in practice. However, I don’t know of an optimized implementation and I would guess that doing a lazy implementation would be slower than computing eigenvalues using LAPACK. As with unpivoted Cholesky, if you have a negative pivot, the matrix is close to an indefinite matrix.