C/C++ Users Journal September, 2005
Cofactor expansion is a direct and well-structured algorithm for computing the determinant of a square matrix of any dimension. Cofactor expansion is defined as the sum of the products of each entry along a row or column of a matrix and the corresponding cofactor of that entry. For an N×N matrix A, the cofactor expansion along the ith row of the matrix is given in Figure 1, where Cij is the cofactor of entry aij, and Mij is the determinant of the submatrix that remains after the ith row and the jth column are removed from the original matrix A.
The algorithm is well structured because it defines a simple recursive pattern. However, the algorithm presents a poor asymptotic runtime. At each stage n in [0,N-1], the expansion sums the product of a subdeterminant of order n-1 (the cofactors) with an entry in the matrix. Thus, the recursive relation is Det(N) = N * Det(N-1) and expands to O(N!).
Despite poor asymptotic performance, cofactor expansion provides a direct method for the computation of the determinant for any N×N matrix. The algorithm can be unrolled to an explicit equation of the matrix coefficients and implemented with no algorithmic overhead. With an efficient way to code the direct computation of the determinant, cofactor expansion would be useful for matrices of small dimensions.
The difficult part of writing this algorithm efficiently is making the computation in-place, and without auxiliary data structures or storage. Each cofactor calculation requires knowledge of a working submatrix. By storing the submatrix indices in compile-time type structures, a metaprogram can drive determinant algorithm code generation for no runtime overhead.
Modern C++ Design [1] presented manipulation of compile-time lists of types with Typelists. Volodymyr Myrnyy extended this idea to Numlists [2], compile-time lists of integers. Numlists are perfect for storing the valid row and column indices of a submatrix, and Numlist type manipulations can be used to access index values and drive template recursion.
The class GenNLSequence (Listing 1) creates a sequential Numlist type structure representing the numbered sequence [0,N-1]. This sequence corresponds to the valid index range of an N×N matrix. The NL::NumAt [2] metafunction provides compile-time access to the index value and the NL::Erase metafunction provides reduction of the index set for template recursion.
The ComputeDet class (Listing 2) is the client interface to the determinant calculation. Its template parameters determine the matrix coefficient type, T, and the dimension of the square input matrix, N. The interface presented is implemented with a vector class specific to a templatized vector library, VectorN<T,N>. However, this interface could be adjusted to your own matrix representations or parameterized through the template. With access to the matrix elements aside, I assume in this article a two-dimensional row major form for its matrices.
ComputeDet provides a single static function, Evaluate(), which takes an array of vectors representing the rows of a matrix. Because cofactor expansion works across rows or columns, the algorithm presented here works if the vectors represent the rows or columns of a matrix. Although the array length could be left undetermined, I prefer specifying the number of vectors expected, designated by the N parameter. The Evaluate() interface is repeated across all auxiliary classes.
To compute the cofactor expansion, ComputeDet forwards the work to the less user friendly, but more powerful, DetEvaluator class.
The DetEvaluator class (Listing 3) is designed to use Numlist type structures to manage submatrix indexing into the vector array. The template parameters MajNListT and MinNListT are Numlist type structures that represent the submatrix indices to be used in the determinant computation. MajNListT manages the first index in two-dimensional matrix addressing, and MinNListT manages the second. Both MajNListT and MinNListT must contain valid indices for an N×N matrix, and must be the same length to represent a square submatrix.
For the full matrix determinant computation, ComputeDet initializes DetEvaluator by passing the full valid index range for an N×N matrix, [0,N-1], generated with GenNLSequence<N>. However, by passing reduced index sequences, DetEvaluator is able to represent submatrix determinant evaluation. The length of the Numlist structures determines the order of the determinant evaluated, and the template parameter N determines the dimension of the matrix interfaced. Because this class can compute any submatrix determinant of a given matrix, it sets up nicely for the recursive evaluation of cofactor expansion.
For cofactor expansion, the DetEvaluator needs to sum the product of the row entries of the submatrix with the corresponding cofactors. DetEvaluator forwards this work to the metaloop class CofactorRecursor.
The CofactorRecursor class (Listing 4) uses a metaloop technique to sum the cofactor products across an entire submatrix row. To do this, CofactorRecursor adds an additional template parameter, MinListIdxN, to track the template recursion. The MinListIdxN parameter is the variable needed to move the summation across a row of the submatrix. It is an index into the MinNListT Numlist.
The major and minor indices are needed to access the matrix entry used in the cofactor product. Because the row to sum across is arbitrary, the first row of the submatrix is always chosen for the major index, MajNListT::value. The minor index is obtained from the MinNListT at the MinListIdxN location, NL::NumAt<MinNListT, MinListIdxN>::value.
To compute the cofactor at each entry, CofactorRecursor reuses the DetEvaluator class, taking advantage of its ability to compute submatrix determinants in-place. To specify the submatrix Numlists for the recursion, CofactorRecursor erases the current major and minor indices from the MajNListT and MinNListT, respectively. The recursive use of DetEvaluator is terminated by the specialization for the case that the submatrix specified is of dimension 1.
After the submatrix determinant is calculated, the cofactor is signed according to the indices of its entry and template recursion is used to accumulate the cofactor products, MinListIdxN is decremented, and template recursion continues until its value is 0.
Finally, the return value of the call to CofactorRecursor by the DetEvaluator class is the determinant value of the submatrix specified by the rows and columns listed in the MajNListT and MinNListT template parameters. Listing 5 shows how easy it is to put these classes to use for the computation of a full matrix determinant.
Numlist manipulation and compile-time metaloops are used to drive template recursion for an efficient implementation of cofactor expansion. Numlists provide the compile-time index list structures that allow processing to happen in place, and the compile-time nature of template recursion is key to coding this algorithm with no runtime overhead.
This implementation (available at http://www.cuj.com/code/) is efficient and generic, capable of calculating the determinant of any submatrix of any order from a matrix of any size in-place. While cofactor expansion is not a good algorithm for large matrices, for small matrices, this algorithm can be used for matrix determinant and inverse functions. At the least, this is an interesting way to calculate an interesting algorithm.