Exploring Quadratic forms with Julia

A quadratic form is a type of polynomial in many variables where all the terms have an exponent of two.

Definition of quadratic form : A quadratic form in $n$ variables $x_1,\dots,x_n$ is given by the expression $q(x) = \sum_{i=1}^{i = n}\sum_{j=1}^{i = n}b_{ij}x_ix_j$ . The quadratic form $q(x)$ can also be expressed in terms of matrix vector product as $q(x) = x^T A x$ where $x = [x_1,\dots,x_n]^T$ and $a_{ij} = \frac{b_{ij} + b_{ji}}{2}$ . Obviously $A$ is a symmetric matrix.

Is matrix $A$ unique for the quadratic form? It seems to be yes.

Hypothesis $H$ : Let $r(x) = x^TCx$ be a quadratic form with asymmmetric matrix $C$ . Then $r(x) = x^TBx$ , where $B = \frac{C+C^T}{2}$ The following Julia code tries to take out any exceptions to the hypothesis.

using Symbolics
@variables x y z
for i in 1:10
         C = rand(1:10,3,3)
         B = (A+A')/2
         q = [x y z]*C*[x y z]'
         eq = simplify(q[1],expand=true)
         q1 = [x y z]*B*[x y z]'
         eq1 = simplify(q1[1],expand=true)
         println(isequal(eq,eq1)) # always true
end

Hypothesis $H$ seems to be true.

Let $n$ be the number of variables in a quadratic form $r(x)$ . Then the number of terms $p$ in $q(x)$ is $\frac{n(n+1)}{2}$ . If $p$ is given then $n = (-1 + \sqrt(1+8p))/2)$ .

Julia function to calculate $n$ given $p$ can be f(p) = Int64(trunc((-1 + sqrt(1+8p))/2)).

A quadratic expression $q(x)$ is

1. Positive definite if $q(x) > 0, \forall x \in R^n$ .

2. Positive semi-definite if $q(x) \ge 0, \forall x \in R^n.$

3. Negative definite if $q(x) < 0, \forall x \in R^n$ .

4. Negative semi-definite if $q(x) \le 0, \forall x \in R^n.$

5. Indefinite if $q(x) > 0$ , for some $x \in R^n$ and $q(x) < 0$ , for some $x \in R^n$

Definiteness of a matrix $A$ can be characterized by eigen values $\lambda_1,\dots,\lambda_n$ of $A$ as

• Positive definite if $\lambda_i > 0, \forall i$ .

• Positive semi-definite if $\lambda_i \ge 0, \forall i$ .

• Negative definite if $\lambda_i < 0, \forall i$ .

• Negative semi-definite if $\lambda_i \le 0, \forall i.$

• Indefinite if $\lambda_i > 0$ , for some $i$ and $\lambda_j < 0$ , for some $j, i \neq j$

The following code in Julia detemines the definiteness of a matrix $A$ in terms of eigenvalues.

using LinearAlgebra
function getDefiniteness(A)

    #@assert isequal(A,A') "Matrix is not symmetric"
ishermitian(A) || throw(ArgumentError("Matrix must be Hermitian"))

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

    if isposdef(-A)
        return "Matrix is negative definite"
    end

    #x,y = eigmin(A),eigmax(A)
    z = eigvals(A)
    x,y = z[1],z[end]

    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 x < 0 && y < 0 #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)

By diagonalizing the symmetric matrix associated with a quadratic form, it is possible to eliminate cross terms.

The Principal Axes Theorem is a mathematical concept related to quadratic forms, which states that a quadratic form $q(x)$ given by $q(x) = x^T A x$ can be diagonalized by an orthonormal set of eigenvectors of $A$ , which is called the set of principal axes. The theorem is based on the spectral theorem related to a symmetric matrix, which states that $Q^T A Q = D$ , where $D$ is a diagonal matrix containing the eigenvalues of $A$ in its diagonal and $Q$ is the eigen matrix of $A$ . By substituting $y = Qx$ , the quadratic form can be expressed as $q(x) = y^T D y = \lambda_1 y_1^2 + \cdots + \lambda_n y_n^2$ , where $\lambda_i$ are the eigenvalues of $A$ .

The minimum and maximum of a quadratic form subjected to the constraint $x_1^2 + \cdots + x_n^2 = 1$ can be calculated by principal axes theorem.

Let $q(x) = x^T Ax, ||x|| = 1$ .

Minimum of $q(x) = \lambda_1$ is the minimum of eigenvalues of $A$ .

Maximum of $q(x) = \lambda_n$ is the maximum of eigenvalues of $A$ .

julia> using  LinearAlgebra

julia> A = rand(1:10,4,4)
4×4 Matrix{Int64}:
 5   6  2   7
 9   8  5   4
 2   7  5  10
 9  10  8   2

julia> A = (A+A')/2
4×4 Matrix{Float64}:
 5.0  7.5  2.0  8.0
 7.5  8.0  6.0  7.0
 2.0  6.0  5.0  9.0
 8.0  7.0  9.0  2.0

julia> minim,maxim = eigmin(A),eigmax(A)
(-7.833772750196167, 25.006216466801746)

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

julia> A = Hermitian(rand(100,100));

julia> @btime a,b = eigminmax(A)
  920.200 μs (32 allocations: 145.23 KiB)
(-5.668761039486605, 49.66411909262399)

julia> @btime x,y = eigmin(A),eigmax(A)
  1.602 ms (25 allocations: 229.81 KiB)
(-5.668761039486608, 49.66411909262401)

Discussion here : https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138