# 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:** 1\
**Showing post:** 3

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

---

_[View the full topic](https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138)._