我们知道用solve可以解y %*% x == b的方程
本文要讲的是一种特殊方程的另一种解法, 这种方程的特点是系数矩阵呈现上三角或下三角特征.
这种方程可以使用backsolve&fowardsolve函数来解.
函数参数 :
> args(backsolve)
function (r, x, k = ncol(r), upper.tri = TRUE, transpose = FALSE)
> args(forwardsolve)
function (l, x, k = ncol(l), upper.tri = FALSE, transpose = FALSE)
其中：r或者l为n×n维三角矩阵，x为n×1维向量(计算乘积时自动转换为矩阵)，对给定不同的upper.tri和transpose的值，方程的形式不同 .
这两个函数用法一致, 只是
upper.tri参数的
默认值不一样. 如果它们参数一致, 得到的结果是一致的.
backsolve(r, x, k = ncol(r), upper.tri = TRUE,
transpose = FALSE)
forwardsolve(l, x, k = ncol(l), upper.tri = FALSE,
transpose = FALSE)
对于函数backsolve()而言，
例如：
> A=matrix(1:9,3,3)
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> x=c(1,2,3)
> x
[1] 1 2 3
> B=A
> upper.tri(B)
[,1] [,2] [,3]
[1,] FALSE TRUE TRUE
[2,] FALSE FALSE TRUE
[3,] FALSE FALSE FALSE
> B[upper.tri(B)]=0 # 使用upper.tri将上三角设置为0
> B
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 2 5 0
[3,] 3 6 9
> C=A
> lower.tri(C)
[,1] [,2] [,3]
[1,] FALSE FALSE FALSE
[2,] TRUE FALSE FALSE
[3,] TRUE TRUE FALSE
> C[lower.tri(C)]=0 # 使用lower.tri将下三角设置为0
> C
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 0 5 8
[3,] 0 0 9
> backsolve(A,x,upper.tri=T,transpose=T) # 参数upper.tri=T表示保留上三角, 即下三角为0
[1] 1.00000000 -0.40000000 -0.08888889 # transpose=T 表示对"执行上三角值保留后的值"进行行列转换
结果和solve(t(C),x)一致.
> solve(t(C),x)
[1] 1.00000000 -0.40000000 -0.08888889
> backsolve(A,x,upper.tri=T,transpose=F) # 参数upper.tri=T表示保留上三角, 即下三角为0
[1] -0.8000000 -0.1333333 0.3333333 # transpose=F 表示对"执行上三角值保留后的值"不进行行列转换
结果和solve(C,x)一致.
> solve(C,x)
[1] -0.8000000 -0.1333333 0.3333333
> backsolve(A,x,upper.tri=F,transpose=T) # 参数upper.tri=F表示不保留上三角, 即上三角为0
[1] 1.111307e-17 2.220446e-17 3.333333e-01
> solve(t(B),x)
[1] 1.110223e-17 2.220446e-17 3.333333e-01
> backsolve(A,x,upper.tri=F,transpose=F) # 参数upper.tri=F表示不保留上三角, 即上三角为0
[1] 1 0 0
> solve(B,x)
[1] 1.000000e+00 -1.540744e-33 -1.850372e-17
对于函数forwardsolve()而言，当传入参数与backsolve()一致时, 结果一致.
例如：
[参考]
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> B
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 2 5 0
[3,] 3 6 9
> C
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 0 5 8
[3,] 0 0 9
> x
[1] 1 2 3
> forwardsolve(A,x,upper.tri=T,transpose=T)
[1] 1.00000000 -0.40000000 -0.08888889
> solve(t(C),x)
[1] 1.00000000 -0.40000000 -0.08888889
> forwardsolve(A,x,upper.tri=T,transpose=F)
[1] -0.8000000 -0.1333333 0.3333333
> solve(C,x)
[1] -0.8000000 -0.1333333 0.3333333
> forwardsolve(A,x,upper.tri=F,transpose=T)
[1] 1.111307e-17 2.220446e-17 3.333333e-01
> solve(t(B),x)
[1] 1.110223e-17 2.220446e-17 3.333333e-01
> forwardsolve(A,x,upper.tri=F,transpose=F)
[1] 1 0 0
> solve(B,x)
[1] 1.000000e+00 -1.540744e-33 -1.850372e-17
1. > help(backsolve)
backsolve package:base R Documentation
Solve an Upper or Lower Triangular System
Description:
Solves a triangular system of linear equations.
Usage:
backsolve(r, x, k = ncol(r), upper.tri = TRUE,
transpose = FALSE)
forwardsolve(l, x, k = ncol(l), upper.tri = FALSE,
transpose = FALSE)
Arguments:
r, l: an upper (or lower) triangular matrix giving the coefficients
for the system to be solved. Values below (above) the
diagonal are ignored.
x: a matrix whose columns give the right-hand sides for the
equations.
k: The number of columns of ‘r’ and rows of ‘x’ to use.
upper.tri: logical; if ‘TRUE’ (default), the _upper_ _tri_angular part
of ‘r’ is used. Otherwise, the lower one.
transpose: logical; if ‘TRUE’, solve r' * y = x for y, i.e., ‘t(r) %*%
y == x’.
Details:
Solves a system of linear equations where the coefficient matrix
is upper (or ‘right’, ‘R’) or lower (‘left’, ‘L’) triangular.
‘x <- backsolve (R, b)’ solves R x = b, and
‘x <- forwardsolve(L, b)’ solves L x = b, respectively.
The ‘r’/‘l’ must have at least ‘k’ rows and columns, and ‘x’ must
have at least ‘k’ rows.
This is a wrapper for the level-3 BLAS routine ‘dtrsm’.
Value:
The solution of the triangular system. The result will be a
vector if ‘x’ is a vector and a matrix if ‘x’ is a matrix.
References:
Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) _The New S
Language_. Wadsworth & Brooks/Cole.
Dongarra, J. J., Bunch, J. R., Moler, C. B. and Stewart, G. W.
(1978) _LINPACK Users Guide_. Philadelphia: SIAM Publications.
See Also:
‘chol’, ‘qr’, ‘solve’.
Examples:
## upper triangular matrix 'r':
r <- rbind(c(1,2,3),
c(0,1,1),
c(0,0,2))
( y <- backsolve(r, x <- c(8,4,2)) ) # -1 3 1
r %*% y # == x = (8,4,2)
backsolve(r, x, transpose = TRUE) # 8 -12 -5
