Subsections

# Solving linear equation systems

When dealing with non-linear networks the number of equation systems to be solved depends on the required precision of the solution and the average necessary iterations until the solution is stable. This emphasizes the meaning of the solving procedures choice for different problems.

The equation systems

 (15.58)

solution can be written as

 (15.59)

## Matrix inversion

The elements of the inverse of the matrix are

 (15.60)

whereas is the matrix elements cofactor. The cofactor is the sub determinant (i.e. the minor) of the element multiplied with . The determinant of a square matrix can be recursively computed by either of the following equations.

 using the -th column (15.61) using the -th row (15.62)

This method is called the Laplace expansion. In order to save computing time the row or column with the most zeros in it is used for the expansion expressed in the above equations. A sub determinant -th order of a matrix's element of -th order is the determinant which is computed by cancelling the -th row and -th column. The following example demonstrates calculating the determinant of a 4th order matrix with the elements of the 3rd row.

 (15.63)

This recursive process for computing the inverse of a matrix is most easiest to be implemented but as well the slowest algorithm. It requires approximately operations.

## Gaussian elimination

The Gaussian algorithm for solving a linear equation system is done in two parts: forward elimination and backward substitution. During forward elimination the matrix A is transformed into an upper triangular equivalent matrix. Elementary transformations due to an equation system having the same solutions for the unknowns as the original system.

 (15.64)

The modifications applied to the matrix A in order to achieve this transformations are limited to the following set of operations.

• multiplication of a row with a scalar factor
• addition or subtraction of multiples of rows
• exchanging two rows of a matrix

### Step 1: Forward elimination

The transformation of the matrix A is done in elimination steps. The new matrix elements of the k-th step with are computed with the following recursive formulas.

 and (15.65) and (15.66) (15.67)

The triangulated matrix can be used to calculate the determinant very easily. The determinant of a triangulated matrix is the product of the diagonal elements. If the determinant is non-zero the equation system has a solution. Otherwise the matrix A is singular.

 (15.68)

When using row and/or column pivoting the resulting determinant may differ in its sign and must be multiplied with whereas is the number of row and column substitutions.

### Finding an appropriate pivot element

The Gaussian elimination fails if the pivot element turns to be zero (division by zero). That is why row and/or column pivoting must be used before each elimination step. If a diagonal element , then exchange the pivot row with the row having the coefficient with the largest absolute value. The new pivot row is and the new pivot element is going to be . If no such pivot row can be found the matrix is singular.

Total pivoting looks for the element with the largest absolute value within the matrix and exchanges rows and columns. When exchanging columns in equation systems the unknowns get reordered as well. For the numerical solution of equation systems with Gaussian elimination column pivoting is clever, and total pivoting recommended.

In order to improve numerical stability pivoting should also be applied if because division by small diagonal elements propagates numerical (rounding) errors. This appears especially with poorly conditioned (the two dimensional case: two lines with nearly the same slope) equation systems.

### Step 2: Backward substitution

The solutions in the vector x are obtained by backward substituting into the triangulated matrix. The elements of the solution vector x are computed by the following recursive equations.

 (15.69) (15.70)

The forward elimination in the Gaussian algorithm requires approximately , the backward substitution operations.

## Gauss-Jordan method

The Gauss-Jordan method is a modification of the Gaussian elimination. In each k-th elimination step the elements of the k-th column get zero except the diagonal element which gets 1. When the right hand side vector z is included in each step it contains the solution vector x afterwards.

The following recursive formulas must be applied to get the new matrix elements for the k-th elimination step. The k-th row must be computed first.

 (15.71) (15.72)

Then the other rows can be calculated with the following formulas.

 (15.73) (15.74)

Column pivoting may be necessary in order to avoid division by zero. The solution vector x is not harmed by row substitutions. When the Gauss-Jordan algorithm has been finished the original matrix has been transformed into the identity matrix. If each operation during this process is applied to an identity matrix the resulting matrix is the inverse matrix of the original matrix. This means that the Gauss-Jordan method can be used to compute the inverse of a matrix.

Though this elimination method is easy to implement the number of required operations is larger than within the Gaussian elimination. The Gauss-Jordan method requires approximately operations.

## LU decomposition

LU decomposition (decomposition into a lower and upper triangular matrix) is recommended when dealing with equation systems where the matrix A does not alter but the right hand side (the vector z) does. Both the Gaussian elimination and the Gauss-Jordan method involve both the right hand side and the matrix in their algorithm. Consecutive solutions of an equation system with an altering right hand side can be computed faster with LU decomposition.

The LU decomposition splits a matrix A into a product of a lower triangular matrix L with an upper triangular matrix U.

 with  and (15.75)

The algorithm for solving the linear equation system involves three steps:

• LU decomposition of the coefficient matrix A
• introduction of an (unknown) arbitrary vector y and solving the equation system by forward substitution
• solving the equation system by backward substitution

The decomposition of the matrix A into a lower and upper triangular matrix is not unique. The most important decompositions, based on Gaussian elimination, are the Doolittle, the Crout and the Cholesky decomposition.

If pivoting is necessary during these algorithms they do not decompose the matrix but the product with an arbitrary matrix (a permutation of the matrix ). When exchanging rows and columns the order of the unknowns as represented by the vector changes as well and must be saved during this process for the forward substitution in the algorithms second step.

### Step 1: LU decomposition

Using the decomposition according to Crout the coefficients of the L and U matrices can be stored in place the original matrix A. The upper triangular matrix U has the form

 (15.76)

The diagonal elements are ones and thus the determinant is one as well. The elements of the new coefficient matrix for the k-th elimination step with compute as follows:

 (15.77) (15.78)

Pivoting may be necessary as you are going to divide by the diagonal element .

### Step 2: Forward substitution

The solutions in the arbitrary vector are obtained by forward substituting into the triangulated matrix. At this stage you need to remember the order of unknowns in the vector as changed by pivoting. The elements of the solution vector are computed by the following recursive equation.

 (15.79)

### Step 3: Backward substitution

The solutions in the vector are obtained by backward substituting into the triangulated matrix. The elements of the solution vector are computed by the following recursive equation.

 (15.80)

The division by the diagonal elements of the matrix U is not necessary because of Crouts definition in eq. (15.76) with .

The LU decomposition requires approximately operations for solving a linear equation system. For consecutive solutions the method requires operations.

## QR decomposition

Singular matrices actually having a solution are over- or under-determined. These types of matrices can be handled by three different types of decompositions: Householder, Jacobi (Givens rotation) and singular value decomposition. Householder decomposition factors a matrix into the product of an orthonormal matrix and an upper triangular matrix , such that:

 (15.81)

The Householder decomposition is based on the fact that for any two different vectors, and , with , i.e. different vectors of equal length, a reflection matrix exists such that:

 (15.82)

To obtain the matrix , the vector is defined by:

 (15.83)

The matrix defined by

 (15.84)

is then the required reflection matrix.

The equation system

 (15.85)

With this yields

 (15.86)

Since is triangular the equation system is solved by a simple matrix-vector multiplication on the right hand side and backward substitution.

### Step 1: QR decomposition

Starting with , let = the first column of , and , i.e. a column vector whose first component is the norm of with the remaining components equal to 0. The Householder transformation with will turn the first column of into as with . At each stage , = the kth column of on and below the diagonal with all other components equal to 0, and 's kth component equals the norm of with all other components equal to 0. Letting , the components of the kth column of below the diagonal are each 0. These calculations are listed below for each stage for the matrix A.

 (15.87)

With this first step the upper left diagonal element of the matrix, , has been generated. The elements below are zeroed out. Since can be generated from stored in place of the first column of the multiplication can be performed without actually generating .

 (15.88)

These elimination steps generate the matrix because is orthonormal, i.e.

 (15.89)

After elimination steps the original matrix contains the upper triangular matrix , except for the diagonal elements which can be stored in some vector. The lower triangular matrix contains the Householder vectors .

 (15.90)

With this representation contains both the and matrix, in a packed form, of course: as a composition of Householder vectors and in the upper triangular part and its diagonal vector .

### Step 2: Forming the new right hand side

In order to form the right hand side let remember eq. (15.84) denoting the reflection matrices used to compute .

 (15.91)

Thus it is possible to replace the original right hand side vector by

 (15.92)

which yields for each the following expression:

 (15.93)

The latter is a simple scalar product of two vectors. Performing eq. (15.93) for each Householder vector finally results in the new right hand side vector .

### Step 3: Backward substitution

The solutions in the vector are obtained by backward substituting into the triangulated matrix. The elements of the solution vector are computed by the following recursive equation.

 (15.94)

### Motivation

Though the QR decomposition has an operation count of (which is about six times more than the LU decomposition) it has its advantages. The QR factorization method is said to be unconditional stable and more accurate. Also it can be used to obtain the minimum-norm (or least square) solution of under-determined equation systems.

The circuit in fig. 15.3 has the following MNA representation:

 (15.95)

The second and third row of the matrix are linear dependent and the matrix is singular because its determinant is zero. Depending on the right hand side , the equation system has none or unlimited solutions. This is called an under-determined system. The discussed QR decomposition easily computes a valid solution without reducing accuracy. The LU decomposition would probably fail because of the singularity.

### Least norm problem

With some more effort it is possible to obtain the minimum-norm solution of this problem. The algorithm as described here would probably yield the following solution:

 (15.96)

This is one out of unlimited solutions. The following short description shows how it is possible to obtain the minimum-norm solution. When decomposing the transposed problem

 (15.97)

the minimum-norm solution is obtained by forward substitution of

 (15.98)

and multiplying the result with .

 (15.99)

In the example above this algorithm results in a solution vector with the least vector norm possible:

 (15.100)

This algorithm outline is also sometimes called LQ decomposition because of being a lower triangular matrix used by the forward substitution.

## Singular value decomposition

Very bad conditioned (ratio between largest and smallest eigenvalue) matrices, i.e. nearly singular, or even singular matrices (over- or under-determined equation systems) can be handled by the singular value decomposition (SVD). This type of decomposition is defined by

 (15.101)

where the matrix consists of the orthonormalized eigenvectors associated with the eigenvalues of , consists of the orthonormalized eigenvectors of and is a matrix with the singular values of (non-negative square roots of the eigenvalues of ) on its diagonal and zeros otherwise.

 (15.102)

The singular value decomposition can be used to solve linear equation systems by simple substitutions

 (15.103) (15.104) (15.105)

since

 (15.106)

To obtain the decomposition stated in eq. (15.101) Householder vectors are computed and their transformations are applied from the left-hand side and right-hand side to obtain an upper bidiagonal matrix which has the same singular values as the original matrix because all of the transformations introduced are orthogonal.

 (15.107)

Specifically, annihilates the subdiagonal elements in column and zeros out the appropriate elements in row .

 (15.108)

Afterwards an iterative process (which turns out to be a QR iteration) is used to transform the bidiagonal matrix into a diagonal form by applying successive Givens transformations (therefore orthogonal as well) to the bidiagonal matrix. This iteration is said to have cubic convergence and yields the final singular values of the matrix .

 (15.109)

 (15.110)

Each of the transformations applied to the bidiagonal matrix is also applied to the matrices and which finally yield the and matrices after convergence.

So far for the algorithm outline. Without the very details the following sections briefly describe each part of the singular value decomposition.

### Notation

Beforehand some notation marks are going to be defined.

• Conjugate transposed (or adjoint):

• Euclidean norm:

• Hermitian (or self adjoint):

whereas denotes the conjugate transposed matrix of . In the real case the matrix is then said to be symmetric''.

• Unitary:

Real matrices with this property are called orthogonal''.

### Householder reflector

A Householder matrix is an elementary unitary matrix that is Hermitian. Its fundamental use is their ability to transform a vector to a multiple of , the first column of the identity matrix. The elementary Hermitian (i.e. the Householder matrix) is defined as

 (15.111)

Beside excellent numerical properties, their application demonstrates their efficiency. If is a matrix, then

 (15.112)

and hence explicit formation and storage of is not required. Also columns (or rows) can be transformed individually exploiting the fact that yields a scalar product for single columns or rows.

#### Specific case

In order to reduce a 44 matrix to upper triangular form successive Householder reflectors must be applied.

 (15.113)

In the first step the diagonal element gets replaced and its below elements get annihilated by the multiplication with an appropriate Householder vector, also the remaining right-hand columns get modified.

 (15.114)

This process must be repeated

 (15.115)

 (15.116)

 (15.117)

until the matrix contains an upper triangular matrix . The matrix can be expressed as the the product of the Householder vectors. The performed operations deliver

 (15.118)

since is unitary. The matrix itself can be expressed in terms of using the following transformation.

 (15.119) (15.120) (15.121)

The eqn. (15.119)-(15.121) are necessary to be mentioned only in case is not Hermitian, but still unitary. Otherwise there is no difference computing or using the Householder vectors. No care must be taken in choosing forward or backward accumulation.

#### General case

In the general case it is necessary to find an elementary unitary matrix

 (15.122)

that satisfies the following three conditions.

 (15.123)

When choosing the elements it is possible the store the Householder vectors as well as the upper triangular matrix in the same storage of the matrix . The Householder matrices can be completely restored from the Householder vectors.

 (15.124)

There exist several approaches to meet the conditions expressed in eq. (15.123). For fewer computational effort it may be convenient to choose to be real valued. With the notation

 (15.125)

one possibility is to define the following calculation rules.

 (15.126) (15.127) (15.128) (15.129) (15.130)

These definitions yield a complex , thus is no more Hermitian but still unitary.

 (15.131)

### Givens rotation

A Givens rotation is a plane rotation matrix. Such a plane rotation matrix is an orthogonal matrix that is different from the identity matrix only in four elements.

 (15.132)

The elements are usually chosen so that

 (15.133)

The most common use of such a plane rotation is to choose and such that for a given and

 (15.134)

multiplication annihilates the lower vector entry. In such an application the matrix is often termed Givens rotation'' matrix. The following equations satisfy eq. (15.134) for a given and exploiting the conditions given in eq. (15.133).

 (15.135)

 (15.136)

### Eigenvalues of a 2-by-2 matrix

The eigenvalues of a 2-by-2 matrix

 (15.137)

can be obtained directly from the quadratic formula. The characteristic polynomial is

 (15.138)

The polynomial yields the two eigenvalues.

 (15.139)

For a symmetric matrix (i.e. ) eq.(15.139) can be rewritten to:

 (15.140)

### Step 1: Bidiagonalization

In the first step the original matrix is bidiagonalized by the application of Householder reflections from the left and right hand side. The matrices and can each be determined as a product of Householder matrices.

 (15.141)

Each of the required Householder vectors are created and applied as previously defined. Suppose a matrix, then applying the first Householder vector from the left hand side eliminates the first column and yields

 (15.142) Next, a Householder vector is applied from the right hand side to annihilate the first row. (15.143) Again, a Householder vector is applied from the left hand side to annihilate the second column. (15.144) This process is continued until (15.145)

For each of the Householder transformations from the left and right hand side the appropriate values must be stored in separate vectors.

### Step 2: Matrix reconstructions

Using the Householder vectors stored in place of the original matrix and the appropriate value vectors it is now necessary to unpack the and matrices. The diagonal vector and the super-diagonal vector can be saved in separate vectors previously. Thus the matrix can be unpacked in place of the matrix and the matrix is unpacked in a separate matrix.

There are two possible algorithms for computing the Householder product matrices, i.e. forward accumulation and backward accumulation. Both start with the identity matrix which is successively multiplied by the Householder matrices either from the left or right.

 (15.146) (15.147)

Recall that the leading portion of each Householder matrix is the identity except the first. Thus, at the beginning of backward accumulation, is mostly the identity'' and it gradually becomes full as the iteration progresses. This pattern can be exploited to reduce the number of required flops. In contrast, is full in forward accumulation after the first step. For this reason, backward accumulation is cheaper and the strategy of choice. When unpacking the matrix in place of the original matrix it is necessary to choose backward accumulation anyway.

 (15.148) (15.149)

Unpacking the matrix is done in a similar way also performing successive Householder matrix multiplications using backward accumulation.

### Step 3: Diagonalization - shifted QR iteration

At this stage the matrices and exist in unfactored form. Also there are the diagonal vector and the super-diagonal vector . Both vectors are real valued. Thus the following algorithm can be applied even though solving a complex equation system.

 (15.150)

The remaining problem is thus to compute the SVD of the matrix . This is done applying an implicit-shift QR step to the tridiagonal matrix which is a symmetric. The matrix is not explicitly formed that is why a QR iteration with implicit shifts is applied.

After bidiagonalization we have a bidiagonal matrix :

 (15.151)

The presented method turns into a matrix by applying a set of orthogonal transforms

 (15.152)

The orthogonal matrices and are chosen so that is also a bidiagonal matrix, but with the super-diagonal elements smaller than those of . The eq.(15.152) is repeated until the non-diagonal elements of become smaller than and can be disregarded.

The matrices and are constructed as

 (15.153)

and similarly where and are matrices of simple rotations as given in eq.(15.132). Both and are products of Givens rotations and thus perform orthogonal transforms.

#### Single shifted QR step.

The left multiplication of by replaces two rows of by their linear combinations. The rest of is unaffected. Right multiplication of by similarly changes only two columns of .

A matrix is chosen the way that

 (15.154)

is a QR transform with a shift. Note that multiplying by gives rise to a non-zero element which is below the main diagonal.

 (15.155)

A new rotation angle is then chosen so that multiplication by gets rid of that element. But this will create a non-zero element which is right beside the super-diagonal.

 (15.156)

Then is made to make it disappear, but this leads to another non-zero element below the diagonal, etc.

 (15.157)

In the end, the matrix becomes bidiagonal again. However, because of a special choice of (QR algorithm), its non-diagonal elements are smaller than those of .

Please note that each of the transforms must also be applied to the unfactored and matrices which turns them successively into and

#### Computation of the Wilkinson shift.

For a single QR step the computation of the eigenvalue of the trailing 2-by-2 submatrix of that is closer to the matrix element is required.

 (15.158) (15.159)

The required eigenvalue is called Wilkinson shift, see eq.(15.140) for details. The sign for the eigenvalue is chosen such that it is closer to .

 (15.160) (15.161) (15.162)

whereas

 (15.163)

The Givens rotation is chosen such that

 (15.164)

The special choice of this first rotation in the single QR step ensures that the super-diagonal matrix entries get smaller. Typically, after a few of these QR steps, the super-diagonal entry becomes negligible.

#### Zeros on the diagonal or super-diagonal.

The QR iteration described above claims to hold if the underlying bidiagonal matrix is unreduced, i.e. has no zeros neither on the diagonal nor on the super-diagonal.

When there is a zero along the diagonal, then premultiplication by a sequence of Givens transformations can zero the right-hand super-diagonal entry as well. The inverse rotations must be applied to the matrix.

Thus the problem can be decoupled into two smaller matrices and . The diagonal matrix is successively getting larger for each super-diagonal entry being neglected after the QR iterations.

 (15.165)

Matrix has non-zero super-diagonal entries. If there is any zero diagonal entry in , then the super-diagonal entry can be annihilated as just described. Otherwise the QR iteration algorithm can be applied to .

When there are only matrix entries left (diagonal entries only) the algorithm is finished, then the matrix has been transformed into the singular value matrix .

### Step 4: Solving the equation system

It is straight-forward to solve a given equation system once having the singular value decomposition computed.

 (15.166) (15.167) (15.168) (15.169) (15.170)

The inverse of the diagonal matrix yields

 (15.171)

With being the i-th row of the matrix , the i-th column of the matrix and the i-th singular value eq. (15.170) can be rewritten to

 (15.172)

It must be mentioned that very small singular values corrupt the complete result. Such values indicate (nearly) singular (ill-conditioned) matrices . In such cases, the solution vector obtained by zeroing the small 's and then using equation (15.170) is better than direct-method solutions (such as LU decomposition or Gaussian elimination) and the SVD solution where the small 's are left non-zero. It may seem paradoxical that this can be so, since zeroing a singular value corresponds to throwing away one linear combination of the set of equations that is going to be solved. The resolution of the paradox is that a combination of equations that is so corrupted by roundoff error is thrown away precisely as to be at best useless; usually it is worse than useless since it "pulls" the solution vector way off towards infinity along some direction that is almost a nullspace vector.

## Jacobi method

This method quite simply involves rearranging each equation to make each variable a function of the other variables. Then make an initial guess for each solution and iterate. For this method it is necessary to ensure that all the diagonal matrix elements are non-zero. This is given for the nodal analysis and almostly given for the modified nodal analysis. If the linear equation system is solvable this can always be achieved by rows substitutions.

The algorithm for performing the iteration step writes as follows.

 (15.173)

This has to repeated until the new solution vectors deviation from the previous one is sufficiently small.

The initial guess has no effect on whether the iterative method converges or not, but with a good initial guess (as possibly given in consecutive Newton-Raphson iterations) it converges faster (if it converges). To ensure convergence the condition

 (15.174)

and at least one case

 (15.175)

must apply. If these conditions are not met, the iterative equations may still converge. If these conditions are met the iterative equations will definitely converge.

Another simple approach to a convergence criteria for iterative algorithms is the Schmidt and v. Mises criteria.

 (15.176)

## Gauss-Seidel method

The Gauss-Seidel algorithm is a modification of the Jacobi method. It uses the previously computed values in the solution vector of the same iteration step. That is why this iterative method is expected to converge faster than the Jacobi method.

The slightly modified algorithm for performing the iteration step writes as follows.

 (15.177)

The remarks about the initial guess as well as the convergence criteria noted in the section about the Jacobi method apply to the Gauss-Seidel algorithm as well.

## A comparison

There are direct and iterative methods (algorithms) for solving linear equation systems. Equation systems with large and sparse matrices should rather be solved with iterative methods.

 method precision application programming effort computing complexity notes Laplace expansion numerical errors general straight forward very time consuming Gaussian elimination numerical errors general intermediate Gauss-Jordan numerical errors general intermediate computes the inverse besides LU decomposition numerical errors general intermediate useful for consecutive solutions QR decomposition good general high Singular value decomposition good general very high ill-conditioned matrices can be handled Jacobi very good diagonally dominant systems easy in each iteration step possibly no convergence Gauss-Seidel very good diagonally dominant systems easy in each iteration step possibly no convergence

This document was generated by Stefan Jahn on 2007-12-30 using latex2html.