The Symmetric Indefinite Decomposition
The symmetric indefinite decomposition or Bunch-Kaufman decomposition is defined for symmetric matrices that may not be positive definite. It expresses a matrix as the product of a lower triangular matrix, a block diagonal matrix, and the transpose of the triangular matrix, with pivot permutations.
The symmetric indefinite decomposition uses symmetry to avoid some of the work of general LU factorization. Pivoting allows it to handle symmetric matrices that are not positive definite.
Working with symmetric indefinite decompositions
The symmetric indefinite decomposition is implemented by the SymmetricIndefiniteDecomposition<T> class. For dense symmetric matrices, obtain a decomposition by calling the GetSymmetricIndefiniteDecomposition method on the matrix. This method has two overloads. The first overload has no arguments. The second overload takes a Boolean value that specifies whether the contents of the matrix may be overwritten by the decomposition. The default is false.
var A = Matrix.CopyFromSymmetric(new double[] {
4.16,-3.12, 0.56,-0.10,
0 ,-5.03,-0.83, 1.18,
0 , 0 , 0.76, 0.34,
0 , 0 , 0 , 1.18
}, 4, MatrixTriangle.Upper, MatrixElementOrder.RowMajor);
var ldlt = A.GetSymmetricIndefiniteDecomposition();The Decompose method performs the actual decomposition. The dense implementation copies the matrix if necessary and calls the appropriate LAPACK routine. This method is called by other methods as needed. You will rarely need to call it explicitly.
Once the decomposition is computed, a number of operations can be performed in much less time. You can repeatedly solve a system of linear equations with different right-hand sides. You can also calculate the determinant and the inverse of the base matrix:
var b = Matrix.CopyFrom(new double[] {
8.70,-13.35, 1.89,-4.14,
8.30, 2.13, 1.61, 5.00
}, 4, 2, MatrixElementOrder.ColumnMajor);
var x = ldlt.Solve(b);
var invA = ldlt.GetInverse();
var detA = ldlt.GetDeterminant();The LowerTriangularFactor property returns the lower component, with pivot permutations included. The BlockDiagonalFactor property returns the block diagonal factor. The three entries in Factors give the factors in multiplication order: the first component, the block diagonal component, and the transpose of the first component (or its conjugate transpose for a Hermitian matrix). FactorTriangle identifies whether the first component is upper or lower before pivot permutations.
var L = ldlt.LowerTriangularFactor;
var LT = ldlt.UpperTriangularFactor;Sparse symmetric indefinite matrices
The same GetSymmetricIndefiniteDecomposition factory handles dense and logically symmetric sparse matrices. For a real double-precision sparse compressed column matrix, it returns a sparse implementation through the shared SymmetricIndefiniteDecomposition<T> contract. The concrete sparse implementation is internal. Other numeric types keep the existing dense implementation. This supports symmetric matrices with positive, negative, or zero eigenvalues. Numerical pivoting chooses scalar or two-by-two blocks; a fill-reducing ordering is enabled.
The common factory requires the logical matrix to be symmetric, with matching entries stored on both sides of the diagonal. The matrix must be square and all entries must be finite. Factorization is lazy: it runs on the first solve, factor access, or numerical diagnostic.
SparseCompressedColumnMatrix<double> matrix = Matrix.CreateSparse<double>(3, 3);
// Store both sides for the common matrix factory.
matrix[0, 0] = 4;
matrix[1, 0] = matrix[0, 1] = 1;
matrix[1, 1] = -2;
matrix[2, 1] = matrix[1, 2] = 1;
matrix[2, 2] = 3;
SymmetricIndefiniteDecomposition<double> factor =
matrix.GetSymmetricIndefiniteDecomposition();
Vector<double> solution = factor.Solve(Vector.Create(2.0, 8.0, 7.0));
// solution is [1, -2, 3]. Reuse factor for other right-hand sides.
int positive = factor.PositiveEigenvalueCount; // 2
int negative = factor.NegativeEigenvalueCount; // 1
double condition = factor.EstimateConditionNumber();The matrix is always preserved. Changes made before the first factorization are observed; later changes do not refresh the factors. Create a new decomposition for a changed matrix. Reuse a completed decomposition for multiple right-hand sides, but synchronize access if sharing it across threads.
The three entries in Factors are the first factor, the block diagonal factor, and the transpose of the first factor. Pivot permutations are absorbed into the first factor, so it is generally not triangular in the original coordinates. Their product is the full symmetric matrix, including reflected entries.
Dense and sparse decompositions share the positive, negative, and zero eigenvalue counts and the two-by-two pivot count. These describe the computed pivot blocks without a near-zero tolerance. A singular decomposition can report factors and rank, but solving or inverting it throws a MatrixSingularException. Rank uses a tolerance on pivot-block eigenvalues; use a singular value decomposition when the tolerance must apply to the original matrix's singular values.
Condition-number estimation uses repeated factor solves and linear workspace, without creating an inverse. Requesting the inverse explicitly produces a dense matrix. Dynamic sparse fill can consume substantial memory; this implementation does not guarantee that factors remain sparse when the input is sparse.