15 Linear algebra tools used in the index reduction process
In the index reduction process, to separate the differential and algebraic components, we need to compute the kernel (null space) of \bm{E}, as it represents the algebraic constraints. There are different approaches available to compute it:
LU decomposition/factorization: a decomposition that factors a matrix A as the product of two matrices, \bm{A} = \bm{L}\bm{U}, where \bm{L} is a lower triangular matrix and \bm{U} is an upper triangular matrix;
Singular Value Decomposition (SVD): a factorization that factors a matrix as the product of three matrices \bm{A} = \bm{U} \bm{D} \bm{V}^T where the columns of \bm{U} and \bm{V} are orthonormal and the matrix \bm{D} is diagonal with positive real entries, to analyze matrix rank and separate the space into differential and algebraic components.
15.1 Index Reduction through Kernel Computation
When the matrix \bm{E} is singular, there exists a nonzero vector \bm{v} \neq \bm{0} such that
\bm{E} \bm{v} = \bm{0}.
All such vectors form a linear subspace called the kernel of \bm{E}, denoted as \operatorname{Kernel}(\bm{E}). This subspace can be spanned by a set of basis vectors \{\bm{v}_i\}. Similarly, there exists a nonzero vector \bm{w} \neq \bm{0} such that
\bm{w}^T \bm{E} = \bm{0}^T,
which defines the kernel of the transpose of \bm{E}, denoted \operatorname{Kernel}(\bm{E}^T). This subspace is spanned by basis vectors \{\bm{w}_i\}.
We focus on the kernel of the transpose, cokernel,
\operatorname{Kernel}(\bm{E}^T) = \{\bm{w} \mid \bm{w}^T \bm{E} = \bm{0}\} = \operatorname{Span}\{\bm{w}_1, \bm{w}_2, \ldots, \bm{w}_p\} = \bm{K}.
Denoting \bm{K} = (\bm{w}_1, \bm{w}_2, \ldots, \bm{w}_p), it follows that
\bm{K}^T \bm{E} = \bm{0}.
If we can compute this kernel matrix \bm{K}(t,\bm{y}) for the time- and state-dependent matrix \bm{E}(t, \bm{y}), we can identify an invariant condition for the DAE.
Consider the DAE system , and define
\bm{K}(t,\bm{y}) = \operatorname{Kernel}(\bm{E}(t,\bm{y})^T) = \operatorname{Cokernel}(\bm{E}).
Multiplying the original DAE by \bm{K}(t,\bm{y})^T, we get:
\bm{K}(t,\bm{y})^T \bm{E}(t,\bm{y}) \bm{y}^\prime = \bm{K}(t,\bm{y})^T \bm{G}(t,\bm{y}).
Since \bm{K}^T \bm{E} = \bm{0}, the left side vanishes, yielding
\bm{0} = \bm{K}(t,\bm{y})^T \bm{G}(t,\bm{y}).
This relation defines an invariant that the solution must satisfy.
Next, we extend the kernel matrix \bm{K}(t,\bm{y}) by finding a complementary matrix \bm{N}(t,\bm{y}) such that the combined matrix
\bm{M}(t,\bm{y}) = \begin{bmatrix} \bm{K}(t,\bm{y}) & \bm{N}(t,\bm{y}) \end{bmatrix}
is nonsingular (this process is often called the completion of the kernel).
Multiplying the original DAE by \bm{M}(t,\bm{y})^T separates the system into algebraic and differential parts:
\bm{M}(t,\bm{y})^T \bm{E}(t,\bm{y}) \bm{y}^\prime = \bm{M}(t,\bm{y})^T \bm{G}(t,\bm{y}),
which can be rewritten as
\begin{bmatrix} \bm{K}(t,\bm{y})^T \\ \bm{N}(t,\bm{y})^T \end{bmatrix} \bm{E}(t,\bm{y}) \bm{y}^\prime = \begin{bmatrix} \bm{K}(t,\bm{y})^T \\ \bm{N}(t,\bm{y})^T \end{bmatrix} \bm{G}(t,\bm{y}).
Manipulating the left side using \bm{K}^T \bm{E} = \bm{0} gives
\begin{bmatrix} \bm{0} \\ \bm{N}(t,\bm{y})^T \bm{E}(t,\bm{y}) \end{bmatrix} \bm{y}^\prime = \begin{bmatrix} \bm{K}(t,\bm{y})^T \bm{G}(t,\bm{y}) \\ \bm{N}(t,\bm{y})^T \bm{G}(t,\bm{y}) \end{bmatrix}.
Define
\bm{\tilde{G}}_1 = \bm{N}(t,\bm{y})^T \bm{G}(t,\bm{y}), \quad \bm{\tilde{G}}_2 = \bm{K}(t,\bm{y})^T \bm{G}(t,\bm{y}).
Differentiating \bm{\tilde{G}}_2 with respect to time,
\frac{d}{dt} \bm{\tilde{G}}_2(t,\bm{y}) = \frac{\partial \bm{\tilde{G}}_2(t,\bm{y})}{\partial t} + \frac{\partial \bm{\tilde{G}}_2(t,\bm{y})}{\partial \bm{y}} \bm{y}^\prime.
The reduced DAE system becomes
\begin{bmatrix} \bm{M}_1(t,\bm{y}) \\ - \frac{\partial \bm{\tilde{G}}_2(t,\bm{y})}{\partial \bm{y}} \end{bmatrix} \bm{y}^\prime = \begin{bmatrix} \bm{\tilde{G}}_1(t,\bm{y}) \\ \frac{\partial \bm{\tilde{G}}_2(t,\bm{y})}{\partial t} \end{bmatrix},
where \bm{M}_1(t,\bm{y}) = \bm{N}(t,\bm{y})^T \bm{E}(t,\bm{y}).
Thus, the DAE can be expressed as the system
\begin{cases} \bm{E}_2(t,\bm{y}) \bm{y}^\prime = \bm{G}_2(t,\bm{y}), \\ \bm{\tilde{G}}_2(t,\bm{y}) = \bm{0}, \end{cases}
with
\bm{E}_2(t,\bm{y}) = \begin{bmatrix} \bm{M}_1(t,\bm{y}) \\ - \frac{\partial \bm{\tilde{G}}_2(t,\bm{y})}{\partial \bm{y}} \end{bmatrix}, \quad \bm{G}_2(t,\bm{y}) = \begin{bmatrix} \bm{\tilde{G}}_1(t,\bm{y}) \\ \frac{\partial \bm{\tilde{G}}_2(t,\bm{y})}{\partial t} \end{bmatrix}.
This reduction process decreases the index of the DAE by one.
In summary, the algorithm proceeds as follows:
- Compute the kernel of \bm{E}^T.
- Complete the kernel to obtain a nonsingular matrix \bm{N}.
- Separate the algebraic and differential parts.
- Differentiate the algebraic constraints.
- Repeat as needed to reduce the index.
Kernel computation methods
The goal is to find a matrix \mathbf{K} whose columns form a basis of \mathrm{Kernel}(\mathbf{E}^T), i.e., they satisfy
\mathbf{K}^T \mathbf{E} = \mathbf{0},
and such that the augmented matrix (\mathbf{K}, \mathbf{N}) is nonsingular (invertible), where \mathbf{N} complements \mathbf{K} to form a full-rank basis of the space.
Two classical methods to compute \mathbf{K} are described below.
Kernel computation using LU decomposition
The LU decomposition factorizes the matrix \mathbf{E} into a product of a lower triangular matrix \mathbf{L} and an upper trapezoidal matrix \mathbf{U}, up to row and column permutations given by permutation matrices \mathbf{P} and \mathbf{Q}:
\mathbf{P} \mathbf{E} \mathbf{Q} = \mathbf{L} \mathbf{U}.
Here, \mathbf{P} and \mathbf{Q} are chosen to ensure numerical stability and reveal the rank structure of \mathbf{E}.
Rearranging, the original matrix \mathbf{E} can be expressed as
\mathbf{E} = \mathbf{P}^T \mathbf{L} \mathbf{U} \mathbf{Q}^T.
Define
\mathbf{M} = \mathbf{U} \mathbf{Q}^T,
so that
\mathbf{E} = \mathbf{P}^T \mathbf{L} \mathbf{M}.
Because \mathbf{U} is upper trapezoidal, \mathbf{M} has a block structure of the form
\mathbf{M} = \begin{bmatrix} \mathbf{M}_1 \\ \mathbf{0} \end{bmatrix},
where \mathbf{M}_1 is full rank.
Multiplying both sides by \mathbf{L}^{-1} \mathbf{P}, we obtain
\mathbf{L}^{-1} \mathbf{P} \mathbf{E} = \mathbf{M} = \begin{bmatrix} \mathbf{M}_1 \\ \mathbf{0} \end{bmatrix}.
Now, by pre-multiplying with the block row selectors (\mathbf{I}, \mathbf{0}) and (\mathbf{0}, \mathbf{I}), where \mathbf{I} is the identity matrix and \mathbf{0} the zero matrix, we isolate the nonzero and zero parts respectively:
(\mathbf{I}, \mathbf{0}) \mathbf{L}^{-1} \mathbf{P} \mathbf{E} = \mathbf{M}_1,
(\mathbf{0}, \mathbf{I}) \mathbf{L}^{-1} \mathbf{P} \mathbf{E} = \mathbf{0}.
From this, we define the matrices \mathbf{K} and \mathbf{N} as
(\mathbf{K}, \mathbf{N}) := \mathbf{P}^T \mathbf{L}^{-T} \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ \mathbf{I} & \mathbf{0} \end{bmatrix},
so that
\mathbf{N}^T = (\mathbf{I}, \mathbf{0}) \mathbf{L}^{-1} \mathbf{P}, \quad \text{and} \quad \mathbf{K}^T = (\mathbf{0}, \mathbf{I}) \mathbf{L}^{-1} \mathbf{P}.
By construction, \mathbf{K}^T \mathbf{E} = \mathbf{0}, and (\mathbf{K}, \mathbf{N}) is nonsingular. Thus, the LU decomposition provides a practical way to compute the kernel of \mathbf{E}.
Kernel computation using Singular Value Decomposition (SVD)
Another standard approach is to use the Singular Value Decomposition (SVD), which factorizes \mathbf{E} into orthogonal matrices \mathbf{U}, \mathbf{V}, and a diagonal matrix \mathbf{\Sigma}:
\mathbf{E} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T,
where
\mathbf{U}^T \mathbf{U} = \mathbf{I}, \quad \mathbf{V}^T \mathbf{V} = \mathbf{I}, \quad \text{and} \quad \mathbf{\Sigma} = \mathrm{diag}(\sigma_1, \ldots, \sigma_r, 0, \ldots, 0).
Partitioning the matrices according to the rank r of \mathbf{E}, we write
\mathbf{E} = \left[ \begin{array}{c|c} \mathbf{U}_{11} & \mathbf{U}_{12} \\ \hline \mathbf{U}_{21} & \mathbf{U}_{22} \end{array} \right] \left[ \begin{array}{c|c} \mathbf{\Sigma}_r & \mathbf{0} \\ \hline \mathbf{0} & \mathbf{0} \end{array} \right] \left[ \begin{array}{c|c} \mathbf{V}_{11}^T & \mathbf{V}_{12}^T \\ \hline \mathbf{V}_{21}^T & \mathbf{V}_{22}^T \end{array} \right].
Multiplying both sides by \mathbf{U}^T and then selecting the lower block with (\mathbf{0}, \mathbf{I}), we get
(\mathbf{0}, \mathbf{I}) \mathbf{U}^T \mathbf{E} = (\mathbf{0}, \mathbf{I}) \left[ \begin{array}{c|c} \mathbf{\Sigma}_r & \mathbf{0} \\ \hline \mathbf{0} & \mathbf{0} \end{array} \right] \left[ \begin{array}{c|c} \mathbf{V}_{11}^T & \mathbf{V}_{12}^T \\ \hline \mathbf{V}_{21}^T & \mathbf{V}_{22}^T \end{array} \right] = \mathbf{0}.
This directly shows that the columns of (\mathbf{0}, \mathbf{I}) \mathbf{U}^T form a basis of the kernel of \mathbf{E}.
Defining
(\mathbf{K}, \mathbf{N}) = \mathbf{U} \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ \mathbf{I} & \mathbf{0} \end{bmatrix},
we have
\mathbf{N}^T = (\mathbf{I}, \mathbf{0}) \mathbf{U}^T, \quad \text{and} \quad \mathbf{K}^T = (\mathbf{0}, \mathbf{I}) \mathbf{U}^T,
and thus \mathbf{K}^T \mathbf{E} = \mathbf{0}.
While the SVD method is robust and numerically stable, it is computationally more expensive than the LU approach and therefore less suited for large-scale or real-time applications. Also for symbolic computations, it is preferable to use LU decomposition.