(* :Name: LinearAlgebra`GaussLagrange` *) (* :Title: Gauss-Lagrange Algorithm *) (* :Author: Gabriel Chenevert Universite de Montreal cheneveg@jsp.umontreal.ca *) (* :Summary: This package implements the Gauss-Lagrange algorithm to find the canonical form under congruence of a symmetric matrix associated with a real quadratic form. This allows one to classify all real quadratic forms, and in particular to determine whether a given quadratic form is positive definite or not. The package also implements elementary row and column operations on any matrix. *) (* :Context: LinearAlgebra`GaussLagrange` *) (* :Package Version: 1.1 *) (* :History: Original version by Gabriel Chenevert, April 1999. Revised and expanded by Gabriel Chenevert, July 1999. *) (* :Keywords: Gauss-Lagrange algorithm, real quadratic form, canonical form of a symmetric matrix under congruence, definiteness *) (* :Source: Any elementary linear algebra textbook. *) (* :Documentation: Notebook GaussLagrange.nb *) (* :Mathematica Version: 3.0 *) (* :Limitations: The core function of the package (GaussLagrange) is written using predominantly a procedural style of programming. Using functional or rule-based programming instead would probably make it more efficient. *) BeginPackage["LinearAlgebra`GaussLagrange`"] Unprotect[AssociatedMatrix, Col, GaussLagrange, R, Summary] Remove[AssociatedMatrix, Col, GaussLagrange, R, Summary] (* Messages *) AssociatedMatrix::usage = "AssociatedMatrix[q, x] gives the matrix associated with the real quadratic form q in vector of variables x, i.e. a symmetric matrix m such that x . m . x == q." AssociatedMatrix::quad = "`1` is not a quadratic form in `2`." Col::usage = "Col[{i, j}, m] swaps the ith and jth columns of the matrix m. Col[i, c, m] multiplies by c the ith column of m. Col[i, j, c, m] adds c times the jth column to the ith column of m." GaussLagrange::usage = "GaussLagrange[m] returns a list {d, p} where d is the canonical form under congruence of the real symmetric matrix m and p is an invertible matrix such that p . m . Transpose[p] == d." GaussLagrange::sym = "Matrix `` is not real and symmetric." R::usage = "R[{i, j}, m] swaps the ith and jth rows of matrix m. R[i, c, m] multiplies by c the ith row of m. R[i, j, c, m] adds c times the jth row to the ith row of m." Summary::usage = "Summary[q, x] prints a summary of the properties of the real quadratic form q in vector of variables x. Summary[m], where m is a real symmetric matrix, uses directly the matrix associated to a quadratic form." Begin["`Private`"] (* Definitions *) $delta[i_Integer, j_Integer] := If[ i == j, 1, 0] $rempl[v_, i_] := Table[ v[[j]] -> $delta[i, j], {j, Length[v]} ] $rempl[v_, k_, l_] := Module[ {m}, m = Table[v[[i]] v[[j]], {i, Length[v]}, {j, Length[v]}]; Flatten[ Table[ m[[i, j]] -> If[{i, j} == {k, l} || {i, j} == {l, k}, 1, 0], {i, Length[m]}, {j, Length[m]} ] ] // Union ] AssociatedMatrix[q_, x_List] := Module[ {res}, res = Table[ If[ i == j, Expand[q] /. $rempl[x, i], Expand[q/2] /. $rempl[x, i, j] ], {i, Length[x]}, {j, Length[x]} ]; If[Simplify[x . res . x - q] =!= 0, Message[AssociatedMatrix::quad, q, x]; res = $Failed ]; res /; res =!= $Failed ] Col[i_Integer, c_, m_List?MatrixQ] := Transpose[R[i, c, Transpose[m]]] Col[{i_Integer, j_Integer}, m_List?MatrixQ] := Transpose[R[{i, j}, Transpose[m]]] Col[i_Integer, j_Integer, c_, m_List?MatrixQ] := Transpose[R[i, j, c, Transpose[m]]] GaussLagrange[a_List?MatrixQ] := Module[ {m = a, n = Length[a], p, temp}, If[m =!= Transpose[m] || !FreeQ[m, Complex], Message[GaussLagrange::sym, m]; p = $Failed; Goto[end]; ]; p = IdentityMatrix[n]; For[j = 1, j <= n, j++, k = j + 1; (* Obtain a pivot != 0 *) While[m[[j, j]] == 0 && k <= n, m = R[j, k, 1, Col[j, k, 1, m]]; p = R[j, k, 1, p]; k++ ]; If[m[[j, j]] == 0, Break[]]; (* Cancel the rest of the row & column *) For[i = j + 1, i <= n, i++, temp = m[[i, j]] / m[[j, j]]; m = R[i, j, -temp, Col[i, j, -temp, m]]; p = R[i, j, -temp, p]; ]; (* Normalize *) temp = Sqrt[Abs[m[[j, j]]]]; If[temp != 0, m = R[j, 1/temp, Col[j, 1/temp, m]]; p = R[j, 1/temp, p]; ]; ]; (* Sort diagonal elements *) For[i = 1, i < n, i++, For[j = i + 1, j <= n, j++, temp = (m[[i, i]] == 0 && m[[j, j]] != 0) || (m[[i, i]] == -1 && m[[j, j]] == 1); If[temp, m = R[{i, j}, Col[{i, j}, m]]; p = R[{i, j}, p]; ]; ]; ]; Label[end]; {m, p} /; p =!= $Failed ] R[i_Integer, c_, m_List?MatrixQ] := Table[If[ k == i, c m[[k,l]], m[[k,l]] ], {k, Length[m]}, {l, Length[m[[1]]]} ] R[i_Integer, j_Integer, c_, m_List?MatrixQ] := Table[If[ k == i, m[[k,l]] + c m[[j,l]], m[[k,l]] ], {k, Length[m]}, {l, Length[m[[1]]]} ] R[{i_Integer, j_Integer}, m_List?MatrixQ] := Part[m, Range[Length[m]] /. {i -> j, j -> i}] Summary[q_, x_List] := Summary[AssociatedMatrix[q, x]] Summary[m_List?MatrixQ] := Module[ {l = GaussLagrange[m], diag, n = Length[m], r, p}, If[ Head[l] == GaussLagrange, Return[] ]; diag = Table[ l[[1, k, k]], {k, n} ]; Print["Order: ", n]; r = Plus @@ Abs /@ diag; Print["Rank: ", r]; p = Length[ Select[diag, # == 1 &] ]; Print["Index: ", p]; Print["Signature: ", 2p - r]; Print[ Which[ r == 0, "Identically 0", p == n, "Positive definite", p == r, "Positive semi-definite", r - p == n, "Negative definite", p == 0, "Negative semi-definite", True, "Indefinite" ], "." ] ] Protect[AssociatedMatrix, Col, GaussLagrange, R, Summary] End[] EndPackage[]