Julia Community 🟣

Vinod V
Vinod V

Posted on • Updated on

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 nn variables x1,,xnx_1,\dots,x_n is given by the expression q(x)=i=1i=nj=1i=nbijxixjq(x) = \sum_{i=1}^{i = n}\sum_{j=1}^{i = n}b_{ij}x_ix_j . The quadratic form q(x)q(x) can also be expressed in terms of matrix vector product as q(x)=xTAxq(x) = x^T A x where x=[x1,,xn]Tx = [x_1,\dots,x_n]^T and aij=bij+bji2a_{ij} = \frac{b_{ij} + b_{ji}}{2} . Obviously AA is a symmetric matrix.

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

Hypothesis HH : Let r(x)=xTCxr(x) = x^TCx be a quadratic form with asymmmetric matrix CC . Then r(x)=xTBxr(x) = x^TBx , where B=C+CT2B = \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
Enter fullscreen mode Exit fullscreen mode

Hypothesis HH seems to be true.

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

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

A quadratic expression q(x)q(x) is

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

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

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

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

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

Definiteness of a matrix AA can be characterized by eigen values λ1,,λn\lambda_1,\dots,\lambda_n of AA as

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

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

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

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

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

The following code in Julia detemines the definiteness of a matrix AA 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)
Enter fullscreen mode Exit fullscreen mode

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)q(x) given by q(x)=xTAxq(x) = x^T A x can be diagonalized by an orthonormal set of eigenvectors of AA , which is called the set of principal axes. The theorem is based on the spectral theorem related to a symmetric matrix, which states that QTAQ=DQ^T A Q = D , where DD is a diagonal matrix containing the eigenvalues of AA in its diagonal and QQ is the eigen matrix of AA . By substituting y=Qxy = Qx , the quadratic form can be expressed as q(x)=yTDy=λ1y12++λnyn2q(x) = y^T D y = \lambda_1 y_1^2 + \cdots + \lambda_n y_n^2 , where λi\lambda_i are the eigenvalues of AA .

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

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

Minimum of q(x)=λ1q(x) = \lambda_1 is the minimum of eigenvalues of AA .

Maximum of q(x)=λnq(x) = \lambda_n is the maximum of eigenvalues of AA .

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)

Enter fullscreen mode Exit fullscreen mode

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

Top comments (0)