Lecture 11: $LU$ Algorithm Variants

The basic procedures within $LU$ factorization can be adapted to account for special structure of the system matrix $𝑨$ or to obtain properties associated with $𝑨$.

1.Determinants

• $𝑨\in {ℝ}^{m\times m}$ a square matrix, $\det \left(𝑨\right)\in ℝ$ is the oriented volume enclosed by the column vectors of $𝑨$ (a parallelipiped)

• Geometric interpretation of determinants

• Determinant calculation rules

• Algebraic definition of a determinant

Definition. The determinant of a square matrix $𝑨=\left(\begin{array}{lll}{𝒂}_1& \dots & {𝒂}_m\end{array}\right)\in {ℝ}^{m\times m}$ is a real number

giving the (oriented) volume of the parallelepiped spanned by matrix column vectors.

• $m=2$

  ${\displaystyle 𝑨=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right),\det \left(𝑨\right)=\left|\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right|}$

• $m=3$

  ${\displaystyle 𝑨=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right),\det \left(𝑨\right)=\left|\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right|}$

• Computation of a determinant with $m=2$

  ${\displaystyle \left|\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right|=a_{11}{a}_{22}-a_{12}{a}_{21}}$

• Computation of a determinant with $m=3$

  ${\displaystyle \begin{array}{rcl}\left|\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right|& =& a_{11}{a}_{22}{a}_{33}+a_{21}{a}_{32}{a}_{13}+a_{31}{a}_{12}{a}_{23}\\ & & -a_{13}{a}_{22}{a}_{31}-a_{23}{a}_{32}{a}_{11}-a_{33}{a}_{12}{a}_{21}\end{array}}$

• Where do these determinant computation rules come from? Two viewpoints

  • Geometric viewpoint: determinants express parallelepiped volumes

  • Algebraic viewpoint: determinants are computed from all possible products that can be formed from choosing a factor from each row and each column

• $m=2$

  Figure 1.

• In two dimensions a “parallelepiped” becomes a parallelogram with area given as

  ${\displaystyle \left(\mathrm{Area}\right)=(\mathrm{Length}\mathrm{of}\mathrm{Base})\times (\mathrm{Length}\mathrm{of}\mathrm{Height})}$

• Take ${𝒂}_1$ as the base, with length $b=||{𝒂}_1||$. Vector ${𝒂}_1$ is at angle ${\phi }_1$ to $x_1$-axis, ${𝒂}_2$ is at angle ${\phi }_2$ to $x_2$-axis, and the angle between ${𝒂}_1$, ${𝒂}_2$ is $\theta ={\phi }_2-{\phi }_1$. The height has length

  ${\displaystyle h=||{𝒂}_2||\sin \theta =||{𝒂}_2||\sin .(.{\phi }_2-{\phi }_1)=||{𝒂}_2||(\sin {\phi }_2\cos {\phi }_1-\sin {\phi }_1\cos {\phi }_2)}$

• Use $\cos {\phi }_1=a_{11}/||{𝒂}_1||$, $\sin {\phi }_1=a_{12}/||{𝒂}_1||$, $\cos {\phi }_2=a_{21}/||{𝒂}_2||$, $\sin {\phi }_2=a_{22}/||{𝒂}_2||$

  ${\displaystyle \left(\mathrm{Area}\right)=||{𝒂}_1||||{𝒂}_2||(\sin {\phi }_2\cos {\phi }_1-\sin {\phi }_1\cos {\phi }_2)=a_{11}{a}_{22}-a_{12}{a}_{21}}$

• The geometric interpretation of a determinant as an oriented volume is useful in establishing rules for calculation with determinants:

  • Determinant of matrix with repeated columns is zero (since two edges of the parallelepiped are identical). Example for $m=3$

    ${\displaystyle \Delta =\left|\begin{array}{ccc}a& a& u\\ b& b& v\\ c& c& w\end{array}\right|=abw+bcu+cav-ubc-vca-wab=0}$

    This is more easily seen using the column notation

    ${\displaystyle \Delta =\det \left(\begin{array}{cccc}{𝒂}_1& {𝒂}_1& {𝒂}_3& \dots \end{array}\right)=0}$

  • Determinant of matrix with linearly dependent columns is zero (since one edge lies in the 'hyperplane' formed by all the others)

• Separating sums in a column (similar for rows)

  ${\displaystyle \det \left(\begin{array}{cccc}{𝒂}_1+{𝒃}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)=\det \left(\begin{array}{cccc}{𝒂}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)+\det \left(\begin{array}{cccc}{𝒃}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)}$

  with ${𝒂}_i,{𝒃}_1\in {ℝ}^m$

• Scalar product in a column (similar for rows)

  ${\displaystyle \det \left(\begin{array}{cccc}\alpha {𝒂}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)=\alpha \det \left(\begin{array}{cccc}{𝒂}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)}$

  with $\alpha \in ℝ$

• Linear combinations of columns (similar for rows)

  ${\displaystyle \det \left(\begin{array}{cccc}{𝒂}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)=\det \left(\begin{array}{cccc}{𝒂}_1& \alpha {𝒂}_1+{𝒂}_2& \dots & {𝒂}_m\end{array}\right)}$

  with $\alpha \in ℝ$.

• A determinant of size $m$ can be expressed as a sum of determinants of size $m-1$ by expansion along a row or column

  ${\displaystyle \begin{array}{rcl}\left|\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \dots & a_{1m}\\ a_{21}& a_{22}& a_{23}& \dots & a_{2m}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& a_{m3}& \dots & a_{mm}\end{array}\right|& =& a_{11}\left|\begin{array}{cccc}{a}_{22}& a_{23}& \dots & a_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m2}& a_{m3}& \dots & a_{mm}\end{array}\right|-\\ & & a_{12}\left|\begin{array}{cccc}{a}_{21}& a_{23}& \dots & a_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m3}& \dots & a_{mm}\end{array}\right|+\\ & & a_{13}\left|\begin{array}{ccccc}{a}_{21}& a_{22}& a_{24}& \dots & a_{2m}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& a_{m4}& \dots & a_{mm}\end{array}\right|-\\ & & \dots \\ & & +{(-1)}^{m+1}{a}_{1m}\left|\begin{array}{cccc}{a}_{21}& a_{23}& \dots & a_{2,m-1}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m3}& \dots & a_{m,m-1}\end{array}\right|\end{array}}$

• The formal definition of a determinant

  ${\displaystyle \det A=\sum _{\sigma \in \Sigma }\nu \left(\sigma \right)a_{1i_1}{a}_{2i_2}\dots a_{m{i}_m}}$

  requires $mm!$ operations, a number that rapidly increases with $m$

• A more economical determinant is to use row and column combinations to create zeros and then reduce the size of the determinant, an algorithm reminiscent of Gauss elimination for systems

  Example:

  ${\displaystyle \left|\begin{array}{ccc}1& 2& 3\\ -1& 0& 1\\ -2& -1& 4\end{array}\right|=\left|\begin{array}{ccc}1& 2& 3\\ 0& 2& 4\\ 0& 3& 10\end{array}\right|=\left|\begin{array}{cc}2& 4\\ 3& 10\end{array}\right|=20-12=8}$

  The first equality comes from linear combinations of rows, i.e. row 1 is added to row 2, and row 1 multiplied by 2 is added to row 3. These linear combinations maintain the value of the determinant. The second equality comes from expansion along the first column

1.1.Cross product

• Consider $u,v\in {ℝ}^3$. We've introduced the idea of a scalar product

  ${\displaystyle u\cdot v=u^{T}v=u_1{v}_1+u_2{v}_2+u_3{v}_3}$

  in which from two vectors one obtains a scalar

• We've also introduced the idea of an exterior product

  ${\displaystyle u{v}^T=\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right)\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)=\left(\begin{array}{ccc}{u}_1{v}_1& u_1{v}_2& u_1{v}_3\\ u_2{v}_1& u_2{v}_2& u_2{v}_3\\ u_3{v}_1& u_3{v}_2& u_3{v}_3\end{array}\right)}$

  in which a matrix is obtained from two vectors

• Another product of two vectors is also useful, the cross product, most conveniently expressed in determinant-like form

  ${\displaystyle u\times v=\left|\begin{array}{ccc}{𝒆}_1& {𝒆}_2& {𝒆}_3\\ u_1& u_2& u_3\\ v_1& v_2& v_3\end{array}\right|=(u_2{v}_3-v_2{u}_3){𝒆}_1+(u_3{v}_1-v_3{u}_1){𝒆}_2+(u_1{v}_2-v_1{u}_2){𝒆}_3}$

2.Structured Matrices

The special structure of a matrix can be exploited to obtain more efficient factorizations. Evaluation of the linear combination $𝑨𝒙=x_1{𝒂}_1+\cdots +x_n{𝒂}_n$ requires $mn$ floating point operations (flops) for $𝑨\in {ℂ}^{m\times n}$. Evaluation of $p$ linear combinations $𝑨𝑿$, $𝑿\in {ℂ}^{n\times p}$ requires $mnp$ flops. If it is possible to evaluate $𝑨𝒙$ with fewer operations, the matrix is said to be structured. Examples include:

• Banded matrices $𝑨=\left[a_{ij}\right]$, $a_{ij}=0$ if $i-j>l$ or $j-i>u$, with $l,u$ denoting the lower and upper bandwidths. If $l=u=0$ the matrix is diagonal. If $l=u=b$ the matrix is said to have bandwidth $B=2b+1$, i.e., for $b=1$, the matrix is tridiagonal, and for $b=2$ the matrix is pentadiagonal. Lower triangular matrices have $u=0$, while upper triangular matrices have $l=0$. The $𝑨𝒙$ product requires $(l+u+1)m$ flops.

• Sparse matrices have $r$ non-zero elements per row or $c$ non-zero elements per column. The $𝑨𝒙$ product requires $rm$ or $cn$ flops

• Circulant matrices $𝑨=\left[a_{ij}\right]$ are sqaure and have $a_{ij}=f(i-j)$, a property that can be exploited to compute $𝑨𝒙$ using $𝒪(m\log m)$ operations

• For square, rank-deficient matrices $𝑨\in {ℂ}^{m\times m}$, $\mathrm{rank}\left(𝑨\right)=r$, $𝑨𝒙$ can be evaluated in $𝒪(km)$ flops

• When $𝑨,𝑿$ are symmetric (hence square), $𝑨𝑿$ requires $𝒪(m^3/2)$ flops instead of $m^3$.

3.Cholesky factorization of positive definite hermitian matrices

3.1.Symmetric matrices, hermitian matrices

Special structure of a matrix is typically associated with underlying symmetries of a particular phenomenon. For example, the law of action and reaction in dynamics (Newton's third law) leads to real symmetric matrices, $𝑨\in {ℝ}^{m\times m}$, ${𝑨}^T=𝑨$. Consider a system of $m$ point masses with nearest-neighbor interactions on the real line where the interaction force depends on relative position. Assume that the force exerted by particle $i+1$ on particle $i$ is linear

with $u_i$ denoting displacement from an equilibrium position. The law of action and reaction then states that

If the same force law holds at all positions, then

The force on particle $i$ is given by the sum of forces from neighboring particles $i-1,i+1$

Introducing $𝒇,𝒖\in {ℝ}^m$, and assuming $u_0=u_{m+1}=0$, the above is stated as

with $𝑲=k\mathrm{diag}\left(\left[\begin{array}{lll}1& -2& 1\end{array}\right]\right)$ is a symmetric matrix, $𝑲={𝑲}^T$, a direct consequence of the law of action and reaction. The matrix $𝑲$ is in this case tridiagonal as a consequence of the assumption of nearest-neighbor interactions. Recall that matrices represent linear mappings, hence

with $𝒇\left(𝒖\right)$ the force-displacement linear mapping, Fig. 2, obtaining the same symmetric, tri-diagonal matrix.

Figure 2. Image of ${𝒆}_i$ through mapping representing a linear force is $𝒇\left({𝒆}_i\right)=k{\left[\begin{array}{lllll}\dots & 1& -2& 1& \dots \end{array}\right]}^T$.

This concept can be extended to complex matrices $𝑨\in {ℂ}^{m\times m}$ through ${𝑨}^{\ast }=𝑨$, in which case $𝑨$ is said to be self-adjoint or hermitian. Again, this property is often associated with desired physical properties, such as the requirement of real observable quantitites in quantum mechanics. Diagonal elements of a hermitian matrix must be real, and for any $𝒙,𝒚\in {ℂ}^m$, the computation

implies for $𝒙=𝒚$ that

hence ${𝒙}^{\ast }𝑨𝒙$ is real.

3.2.Positive-definite matrices

The work (i.e., energy stored in the system) done by all the forces in the above point mass system is

and physical considerations state that $𝒲⩾0$. This leads the following definitions.

Definition. A hermitian matrix $𝑨\in {ℂ}^{m\times m}$ is positive definite if for any non-zero $𝒙\in {ℂ}^m,$ ${𝒙}^{\ast }𝑨𝒙>0$.

Definition. A hermitian matrix $𝑨\in {ℂ}^{m\times m}$ is positive semi-definite if for any non-zero $𝒙\in {ℂ}^m,$ ${𝒙}^{\ast }𝑨𝒙⩾0$.

If $𝑨$ is hermitian positive definite, then so is ${𝑿}^{\ast }𝑨𝑿$ for any $𝑿\in {ℂ}^{m\times n}$. Choosing

gives ${𝑨}_n={𝑿}^{\ast }𝑨𝑿$, the $n^{\mathrm{th}}$ principal submatrix of $𝑨$, itself a hermitian positive definite matrix. Choosing $𝑿={𝒆}_j$shows that the $j^{\mathrm{th}}$ diagonal element of $𝑨$ is positive, $a_{jj}={𝒆}_j^T𝑨{𝒆}_j>0$

3.3.Symmetric factorization of positive-definite hermitian matrices

The structure of a hermitian positive definite matrix $𝑨\in {ℂ}^{m\times m}$, can be preserved by modification of $LU$-factorization. The resulting algorithm is known as Cholesky factorization, and its first stage is stated as

whence $𝑪=𝑩-𝒘{𝒘}^{\ast }/a_{11}$. Repeating the stage-1 step

leads to

The resulting Cholesky algorithm is half as expensive as standard $LU$-factorization.