Representation of a matrix as a sum
In the mathematical discipline of numerical linear algebra , a matrix splitting is an expression which represents a given matrix as a sum or difference of matrices. Many iterative methods (for example, for systems of differential equations ) depend upon the direct solution of matrix equations involving matrices more general than tridiagonal matrices . These matrix equations can often be solved directly and efficiently when written as a matrix splitting. The technique was devised by Richard S. Varga in 1960.[ 1]
Regular splittings [ edit ]
We seek to solve the matrix equation
A
x
=
k
,
{\displaystyle \mathbf {A} \mathbf {x} =\mathbf {k} ,}
{\displaystyle \mathbf {A} \mathbf {x} =\mathbf {k} ,} 1
where A is a given n × ばつ n non-singular matrix, and k is a given column vector with n components. We split the matrix A into
A
=
B
−
C
,
{\displaystyle \mathbf {A} =\mathbf {B} -\mathbf {C} ,}
{\displaystyle \mathbf {A} =\mathbf {B} -\mathbf {C} ,} 2
where B and C are n × ばつ n matrices. If, for an arbitrary n × ばつ n matrix M , M has nonnegative entries, we write M ≥ 0 . If M has only positive entries, we write M > 0 . Similarly, if the matrix M 1 − M 2 has nonnegative entries, we write M 1 ≥ M 2 .
Definition: A = B − C is a regular splitting of A if B −1 ≥ 0 and C ≥ 0 .
We assume that matrix equations of the form
B
x
=
g
,
{\displaystyle \mathbf {B} \mathbf {x} =\mathbf {g} ,}
{\displaystyle \mathbf {B} \mathbf {x} =\mathbf {g} ,} 3
where g is a given column vector, can be solved directly for the vector x . If (2 ) represents a regular splitting of A , then the iterative method
B
x
(
m
+
1
)
=
C
x
(
m
)
+
k
,
m
=
0
,
1
,
2
,
…
,
{\displaystyle \mathbf {B} \mathbf {x} ^{(m+1)}=\mathbf {C} \mathbf {x} ^{(m)}+\mathbf {k} ,\quad m=0,1,2,\ldots ,}
{\displaystyle \mathbf {B} \mathbf {x} ^{(m+1)}=\mathbf {C} \mathbf {x} ^{(m)}+\mathbf {k} ,\quad m=0,1,2,\ldots ,} 4
where x (0) is an arbitrary vector, can be carried out. Equivalently, we write (4 ) in the form
x
(
m
+
1
)
=
B
−
1
C
x
(
m
)
+
B
−
1
k
,
m
=
0
,
1
,
2
,
…
{\displaystyle \mathbf {x} ^{(m+1)}=\mathbf {B} ^{-1}\mathbf {C} \mathbf {x} ^{(m)}+\mathbf {B} ^{-1}\mathbf {k} ,\quad m=0,1,2,\ldots }
{\displaystyle \mathbf {x} ^{(m+1)}=\mathbf {B} ^{-1}\mathbf {C} \mathbf {x} ^{(m)}+\mathbf {B} ^{-1}\mathbf {k} ,\quad m=0,1,2,\ldots } 5
The matrix D = B −1 C has nonnegative entries if (2 ) represents a regular splitting of A .[ 2]
It can be shown that if A −1 > 0 , then
ρ
(
D
)
{\displaystyle \rho (\mathbf {D} )}
{\displaystyle \rho (\mathbf {D} )} < 1, where
ρ
(
D
)
{\displaystyle \rho (\mathbf {D} )}
{\displaystyle \rho (\mathbf {D} )} represents the spectral radius of D , and thus D is a convergent matrix . As a consequence, the iterative method (5 ) is necessarily convergent .[ 3] [ 4]
If, in addition, the splitting (2 ) is chosen so that the matrix B is a diagonal matrix (with the diagonal entries all non-zero, since B must be invertible ), then B can be inverted in linear time (see Time complexity ).
Matrix iterative methods [ edit ]
Many iterative methods can be described as a matrix splitting. If the diagonal entries of the matrix A are all nonzero, and we express the matrix A as the matrix sum
A
=
D
−
U
−
L
,
{\displaystyle \mathbf {A} =\mathbf {D} -\mathbf {U} -\mathbf {L} ,}
{\displaystyle \mathbf {A} =\mathbf {D} -\mathbf {U} -\mathbf {L} ,} 6
where D is the diagonal part of A , and U and L are respectively strictly upper and lower triangular n × ばつ n matrices, then we have the following.
The Jacobi method can be represented in matrix form as a splitting
x
(
m
+
1
)
=
D
−
1
(
U
+
L
)
x
(
m
)
+
D
−
1
k
.
{\displaystyle \mathbf {x} ^{(m+1)}=\mathbf {D} ^{-1}(\mathbf {U} +\mathbf {L} )\mathbf {x} ^{(m)}+\mathbf {D} ^{-1}\mathbf {k} .}
{\displaystyle \mathbf {x} ^{(m+1)}=\mathbf {D} ^{-1}(\mathbf {U} +\mathbf {L} )\mathbf {x} ^{(m)}+\mathbf {D} ^{-1}\mathbf {k} .}[ 5] [ 6] 7
The Gauss–Seidel method can be represented in matrix form as a splitting
x
(
m
+
1
)
=
(
D
−
L
)
−
1
U
x
(
m
)
+
(
D
−
L
)
−
1
k
.
{\displaystyle \mathbf {x} ^{(m+1)}=(\mathbf {D} -\mathbf {L} )^{-1}\mathbf {U} \mathbf {x} ^{(m)}+(\mathbf {D} -\mathbf {L} )^{-1}\mathbf {k} .}
{\displaystyle \mathbf {x} ^{(m+1)}=(\mathbf {D} -\mathbf {L} )^{-1}\mathbf {U} \mathbf {x} ^{(m)}+(\mathbf {D} -\mathbf {L} )^{-1}\mathbf {k} .}[ 7] [ 8] 8
The method of successive over-relaxation can be represented in matrix form as a splitting
x
(
m
+
1
)
=
(
D
−
ω
L
)
−
1
[
(
1
−
ω
)
D
+
ω
U
]
x
(
m
)
+
ω
(
D
−
ω
L
)
−
1
k
.
{\displaystyle \mathbf {x} ^{(m+1)}=(\mathbf {D} -\omega \mathbf {L} )^{-1}[(1-\omega )\mathbf {D} +\omega \mathbf {U} ]\mathbf {x} ^{(m)}+\omega (\mathbf {D} -\omega \mathbf {L} )^{-1}\mathbf {k} .}
{\displaystyle \mathbf {x} ^{(m+1)}=(\mathbf {D} -\omega \mathbf {L} )^{-1}[(1-\omega )\mathbf {D} +\omega \mathbf {U} ]\mathbf {x} ^{(m)}+\omega (\mathbf {D} -\omega \mathbf {L} )^{-1}\mathbf {k} .}[ 9] [ 10] 9
In equation (1 ), let
A
=
(
6
−
2
−
3
−
1
4
−
2
−
3
−
1
5
)
,
k
=
(
5
−
12
10
)
.
{\displaystyle \mathbf {A} ={\begin{pmatrix}6&-2&-3\\-1&4&-2\\-3&-1&5\end{pmatrix}},\quad \mathbf {k} ={\begin{pmatrix}5\\-12\10円\end{pmatrix}}.}
{\displaystyle \mathbf {A} ={\begin{pmatrix}6&-2&-3\\-1&4&-2\\-3&-1&5\end{pmatrix}},\quad \mathbf {k} ={\begin{pmatrix}5\\-12\10円\end{pmatrix}}.} 10
Let us apply the splitting (7 ) which is used in the Jacobi method: we split A in such a way that B consists of all of the diagonal elements of A , and C consists of all of the off-diagonal elements of A , negated. (Of course this is not the only useful way to split a matrix into two matrices.) We have
B
=
(
6
0
0
0
4
0
0
0
5
)
,
C
=
(
0
2
3
1
0
2
3
1
0
)
,
{\displaystyle {\begin{aligned}&\mathbf {B} ={\begin{pmatrix}6&0&0\0円&4&0\0円&0&5\end{pmatrix}},\quad \mathbf {C} ={\begin{pmatrix}0&2&3\1円&0&2\3円&1&0\end{pmatrix}},\end{aligned}}}
{\displaystyle {\begin{aligned}&\mathbf {B} ={\begin{pmatrix}6&0&0\0円&4&0\0円&0&5\end{pmatrix}},\quad \mathbf {C} ={\begin{pmatrix}0&2&3\1円&0&2\3円&1&0\end{pmatrix}},\end{aligned}}} 11
A
−
1
=
1
47
(
18
13
16
11
21
15
13
12
22
)
,
B
−
1
=
(
1
6
0
0
0
1
4
0
0
0
1
5
)
,
{\displaystyle {\begin{aligned}&\mathbf {A^{-1}} ={\frac {1}{47}}{\begin{pmatrix}18&13&16\11円&21&15\13円&12&22\end{pmatrix}},\quad \mathbf {B^{-1}} ={\begin{pmatrix}{\frac {1}{6}}&0&0\\[4pt]0&{\frac {1}{4}}&0\\[4pt]0&0&{\frac {1}{5}}\end{pmatrix}},\end{aligned}}}
{\displaystyle {\begin{aligned}&\mathbf {A^{-1}} ={\frac {1}{47}}{\begin{pmatrix}18&13&16\11円&21&15\13円&12&22\end{pmatrix}},\quad \mathbf {B^{-1}} ={\begin{pmatrix}{\frac {1}{6}}&0&0\\[4pt]0&{\frac {1}{4}}&0\\[4pt]0&0&{\frac {1}{5}}\end{pmatrix}},\end{aligned}}}
D
=
B
−
1
C
=
(
0
1
3
1
2
1
4
0
1
2
3
5
1
5
0
)
,
B
−
1
k
=
(
5
6
−
3
2
)
.
{\displaystyle {\begin{aligned}\mathbf {D} =\mathbf {B^{-1}C} ={\begin{pmatrix}0&{\frac {1}{3}}&{\frac {1}{2}}\\[4pt]{\frac {1}{4}}&0&{\frac {1}{2}}\\[4pt]{\frac {3}{5}}&{\frac {1}{5}}&0\end{pmatrix}},\quad \mathbf {B^{-1}k} ={\begin{pmatrix}{\frac {5}{6}}\\[4pt]-3\\[4pt]2\end{pmatrix}}.\end{aligned}}}
{\displaystyle {\begin{aligned}\mathbf {D} =\mathbf {B^{-1}C} ={\begin{pmatrix}0&{\frac {1}{3}}&{\frac {1}{2}}\\[4pt]{\frac {1}{4}}&0&{\frac {1}{2}}\\[4pt]{\frac {3}{5}}&{\frac {1}{5}}&0\end{pmatrix}},\quad \mathbf {B^{-1}k} ={\begin{pmatrix}{\frac {5}{6}}\\[4pt]-3\\[4pt]2\end{pmatrix}}.\end{aligned}}}
Since B −1 ≥ 0 and C ≥ 0 , the splitting (11 ) is a regular splitting. Since A −1 > 0 , the spectral radius
ρ
(
D
)
{\displaystyle \rho (\mathbf {D} )}
{\displaystyle \rho (\mathbf {D} )} < 1. (The approximate eigenvalues of D are
λ
i
≈
−
0.4599820
,
−
0.3397859
,
0.7997679.
{\displaystyle \lambda _{i}\approx -0.4599820,-0.3397859,0.7997679.}
{\displaystyle \lambda _{i}\approx -0.4599820,-0.3397859,0.7997679.} ) Hence, the matrix D is convergent and the method (5 ) necessarily converges for the problem (10 ). Note that the diagonal elements of A are all greater than zero, the off-diagonal elements of A are all less than zero and A is strictly diagonally dominant .[ 11]
The method (5 ) applied to the problem (10 ) then takes the form
x
(
m
+
1
)
=
(
0
1
3
1
2
1
4
0
1
2
3
5
1
5
0
)
x
(
m
)
+
(
5
6
−
3
2
)
,
m
=
0
,
1
,
2
,
…
{\displaystyle \mathbf {x} ^{(m+1)}={\begin{pmatrix}0&{\frac {1}{3}}&{\frac {1}{2}}\\[4pt]{\frac {1}{4}}&0&{\frac {1}{2}}\\[4pt]{\frac {3}{5}}&{\frac {1}{5}}&0\end{pmatrix}}\mathbf {x} ^{(m)}+{\begin{pmatrix}{\frac {5}{6}}\\[4pt]-3\\[4pt]2\end{pmatrix}},\quad m=0,1,2,\ldots }
{\displaystyle \mathbf {x} ^{(m+1)}={\begin{pmatrix}0&{\frac {1}{3}}&{\frac {1}{2}}\\[4pt]{\frac {1}{4}}&0&{\frac {1}{2}}\\[4pt]{\frac {3}{5}}&{\frac {1}{5}}&0\end{pmatrix}}\mathbf {x} ^{(m)}+{\begin{pmatrix}{\frac {5}{6}}\\[4pt]-3\\[4pt]2\end{pmatrix}},\quad m=0,1,2,\ldots } 12
The exact solution to equation (12 ) is
x
=
(
2
−
1
3
)
.
{\displaystyle \mathbf {x} ={\begin{pmatrix}2\\-1\3円\end{pmatrix}}.}
{\displaystyle \mathbf {x} ={\begin{pmatrix}2\\-1\3円\end{pmatrix}}.} 13
The first few iterates for equation (12 ) are listed in the table below, beginning with x (0) = (0.0, 0.0, 0.0)T . From the table one can see that the method is evidently converging to the solution (13 ), albeit rather slowly.
x
1
(
m
)
{\displaystyle x_{1}^{(m)}}
{\displaystyle x_{1}^{(m)}}
x
2
(
m
)
{\displaystyle x_{2}^{(m)}}
{\displaystyle x_{2}^{(m)}}
x
3
(
m
)
{\displaystyle x_{3}^{(m)}}
{\displaystyle x_{3}^{(m)}}
0.0
0.0
0.0
0.83333
-3.0000
2.0000
0.83333
-1.7917
1.9000
1.1861
-1.8417
2.1417
1.2903
-1.6326
2.3433
1.4608
-1.5058
2.4477
1.5553
-1.4110
2.5753
1.6507
-1.3235
2.6510
1.7177
-1.2618
2.7257
1.7756
-1.2077
2.7783
1.8199
-1.1670
2.8238
As stated above, the Jacobi method (7 ) is the same as the specific regular splitting (11 ) demonstrated above.
Gauss–Seidel method [ edit ]
Since the diagonal entries of the matrix A in problem (10 ) are all nonzero, we can express the matrix A as the splitting (6 ), where
D
=
(
6
0
0
0
4
0
0
0
5
)
,
U
=
(
0
2
3
0
0
2
0
0
0
)
,
L
=
(
0
0
0
1
0
0
3
1
0
)
.
{\displaystyle \mathbf {D} ={\begin{pmatrix}6&0&0\0円&4&0\0円&0&5\end{pmatrix}},\quad \mathbf {U} ={\begin{pmatrix}0&2&3\0円&0&2\0円&0&0\end{pmatrix}},\quad \mathbf {L} ={\begin{pmatrix}0&0&0\1円&0&0\3円&1&0\end{pmatrix}}.}
{\displaystyle \mathbf {D} ={\begin{pmatrix}6&0&0\0円&4&0\0円&0&5\end{pmatrix}},\quad \mathbf {U} ={\begin{pmatrix}0&2&3\0円&0&2\0円&0&0\end{pmatrix}},\quad \mathbf {L} ={\begin{pmatrix}0&0&0\1円&0&0\3円&1&0\end{pmatrix}}.} 14
We then have
(
D
−
L
)
−
1
=
1
120
(
20
0
0
5
30
0
13
6
24
)
,
{\displaystyle {\begin{aligned}&\mathbf {(D-L)^{-1}} ={\frac {1}{120}}{\begin{pmatrix}20&0&0\5円&30&0\13円&6&24\end{pmatrix}},\end{aligned}}}
{\displaystyle {\begin{aligned}&\mathbf {(D-L)^{-1}} ={\frac {1}{120}}{\begin{pmatrix}20&0&0\5円&30&0\13円&6&24\end{pmatrix}},\end{aligned}}}
(
D
−
L
)
−
1
U
=
1
120
(
0
40
60
0
10
75
0
26
51
)
,
(
D
−
L
)
−
1
k
=
1
120
(
100
−
335
233
)
.
{\displaystyle {\begin{aligned}&\mathbf {(D-L)^{-1}U} ={\frac {1}{120}}{\begin{pmatrix}0&40&60\0円&10&75\0円&26&51\end{pmatrix}},\quad \mathbf {(D-L)^{-1}k} ={\frac {1}{120}}{\begin{pmatrix}100\\-335\233円\end{pmatrix}}.\end{aligned}}}
{\displaystyle {\begin{aligned}&\mathbf {(D-L)^{-1}U} ={\frac {1}{120}}{\begin{pmatrix}0&40&60\0円&10&75\0円&26&51\end{pmatrix}},\quad \mathbf {(D-L)^{-1}k} ={\frac {1}{120}}{\begin{pmatrix}100\\-335\233円\end{pmatrix}}.\end{aligned}}}
The Gauss–Seidel method (8 ) applied to the problem (10 ) takes the form
x
(
m
+
1
)
=
1
120
(
0
40
60
0
10
75
0
26
51
)
x
(
m
)
+
1
120
(
100
−
335
233
)
,
m
=
0
,
1
,
2
,
…
{\displaystyle \mathbf {x} ^{(m+1)}={\frac {1}{120}}{\begin{pmatrix}0&40&60\0円&10&75\0円&26&51\end{pmatrix}}\mathbf {x} ^{(m)}+{\frac {1}{120}}{\begin{pmatrix}100\\-335\233円\end{pmatrix}},\quad m=0,1,2,\ldots }
{\displaystyle \mathbf {x} ^{(m+1)}={\frac {1}{120}}{\begin{pmatrix}0&40&60\0円&10&75\0円&26&51\end{pmatrix}}\mathbf {x} ^{(m)}+{\frac {1}{120}}{\begin{pmatrix}100\\-335\233円\end{pmatrix}},\quad m=0,1,2,\ldots } 15
The first few iterates for equation (15 ) are listed in the table below, beginning with x (0) = (0.0, 0.0, 0.0)T . From the table one can see that the method is evidently converging to the solution (13 ), somewhat faster than the Jacobi method described above.
x
1
(
m
)
{\displaystyle x_{1}^{(m)}}
{\displaystyle x_{1}^{(m)}}
x
2
(
m
)
{\displaystyle x_{2}^{(m)}}
{\displaystyle x_{2}^{(m)}}
x
3
(
m
)
{\displaystyle x_{3}^{(m)}}
{\displaystyle x_{3}^{(m)}}
0.0
0.0
0.0
0.8333
-2.7917
1.9417
0.8736
-1.8107
2.1620
1.3108
-1.5913
2.4682
1.5370
-1.3817
2.6459
1.6957
-1.2531
2.7668
1.7990
-1.1668
2.8461
1.8675
-1.1101
2.8985
1.9126
-1.0726
2.9330
1.9423
-1.0479
2.9558
1.9619
-1.0316
2.9708
Successive over-relaxation method [ edit ]
Let ω = 1.1. Using the splitting (14 ) of the matrix A in problem (10 ) for the successive over-relaxation method, we have
(
D
−
ω
L
)
−
1
=
1
12
(
2
0
0
0.55
3
0
1.441
0.66
2.4
)
,
{\displaystyle {\begin{aligned}&\mathbf {(D-\omega L)^{-1}} ={\frac {1}{12}}{\begin{pmatrix}2&0&0\0円.55&3&0\1円.441&0.66&2.4\end{pmatrix}},\end{aligned}}}
{\displaystyle {\begin{aligned}&\mathbf {(D-\omega L)^{-1}} ={\frac {1}{12}}{\begin{pmatrix}2&0&0\0円.55&3&0\1円.441&0.66&2.4\end{pmatrix}},\end{aligned}}}
(
D
−
ω
L
)
−
1
[
(
1
−
ω
)
D
+
ω
U
]
=
1
12
(
−
1.2
4.4
6.6
−
0.33
0.01
8.415
−
0.8646
2.9062
5.0073
)
,
{\displaystyle {\begin{aligned}&\mathbf {(D-\omega L)^{-1}[(1-\omega )D+\omega U]} ={\frac {1}{12}}{\begin{pmatrix}-1.2&4.4&6.6\\-0.33&0.01&8.415\\-0.8646&2.9062&5.0073\end{pmatrix}},\end{aligned}}}
{\displaystyle {\begin{aligned}&\mathbf {(D-\omega L)^{-1}[(1-\omega )D+\omega U]} ={\frac {1}{12}}{\begin{pmatrix}-1.2&4.4&6.6\\-0.33&0.01&8.415\\-0.8646&2.9062&5.0073\end{pmatrix}},\end{aligned}}}
ω
(
D
−
ω
L
)
−
1
k
=
1
12
(
11
−
36.575
25.6135
)
.
{\displaystyle {\begin{aligned}&\mathbf {\omega (D-\omega L)^{-1}k} ={\frac {1}{12}}{\begin{pmatrix}11\\-36.575\25円.6135\end{pmatrix}}.\end{aligned}}}
{\displaystyle {\begin{aligned}&\mathbf {\omega (D-\omega L)^{-1}k} ={\frac {1}{12}}{\begin{pmatrix}11\\-36.575\25円.6135\end{pmatrix}}.\end{aligned}}}
The successive over-relaxation method (9 ) applied to the problem (10 ) takes the form
x
(
m
+
1
)
=
1
12
(
−
1.2
4.4
6.6
−
0.33
0.01
8.415
−
0.8646
2.9062
5.0073
)
x
(
m
)
+
1
12
(
11
−
36.575
25.6135
)
,
m
=
0
,
1
,
2
,
…
{\displaystyle \mathbf {x} ^{(m+1)}={\frac {1}{12}}{\begin{pmatrix}-1.2&4.4&6.6\\-0.33&0.01&8.415\\-0.8646&2.9062&5.0073\end{pmatrix}}\mathbf {x} ^{(m)}+{\frac {1}{12}}{\begin{pmatrix}11\\-36.575\25円.6135\end{pmatrix}},\quad m=0,1,2,\ldots }
{\displaystyle \mathbf {x} ^{(m+1)}={\frac {1}{12}}{\begin{pmatrix}-1.2&4.4&6.6\\-0.33&0.01&8.415\\-0.8646&2.9062&5.0073\end{pmatrix}}\mathbf {x} ^{(m)}+{\frac {1}{12}}{\begin{pmatrix}11\\-36.575\25円.6135\end{pmatrix}},\quad m=0,1,2,\ldots } 16
The first few iterates for equation (16 ) are listed in the table below, beginning with x (0) = (0.0, 0.0, 0.0)T . From the table one can see that the method is evidently converging to the solution (13 ), slightly faster than the Gauss–Seidel method described above.
x
1
(
m
)
{\displaystyle x_{1}^{(m)}}
{\displaystyle x_{1}^{(m)}}
x
2
(
m
)
{\displaystyle x_{2}^{(m)}}
{\displaystyle x_{2}^{(m)}}
x
3
(
m
)
{\displaystyle x_{3}^{(m)}}
{\displaystyle x_{3}^{(m)}}
0.0
0.0
0.0
0.9167
-3.0479
2.1345
0.8814
-1.5788
2.2209
1.4711
-1.5161
2.6153
1.6521
-1.2557
2.7526
1.8050
-1.1641
2.8599
1.8823
-1.0930
2.9158
1.9314
-1.0559
2.9508
1.9593
-1.0327
2.9709
1.9761
-1.0185
2.9829
1.9862
-1.0113
2.9901
^ Varga (1960)
^ Varga (1960 , pp. 121–122)
^ Varga (1960 , pp. 122–123)
^ Varga (1962 , p. 89)
^ Burden & Faires (1993 , p. 408)
^ Varga (1962 , p. 88)
^ Burden & Faires (1993 , p. 411)
^ Varga (1962 , p. 88)
^ Burden & Faires (1993 , p. 416)
^ Varga (1962 , p. 88)
^ Burden & Faires (1993 , p. 371)
Burden, Richard L.; Faires, J. Douglas (1993), Numerical Analysis (5th ed.), Boston: Prindle, Weber and Schmidt , ISBN 0-534-93219-3 .
Varga, Richard S. (1960). "Factorization and Normalized Iterative Methods". In Langer, Rudolph E. (ed.). Boundary Problems in Differential Equations . Madison: University of Wisconsin Press . pp. 121– 142. LCCN 60-60003 .
Varga, Richard S. (1962), Matrix Iterative Analysis , New Jersey: Prentice-Hall , Bibcode :1962mia..book.....V , LCCN 62-21277 .
Key concepts Problems Hardware Software