LU 分解

# LU 分解

(1)
\begin{align} \Bbb { A = P L U } \end{align}

# 範例一

> library(Matrix)

Attaching package: 'Matrix'

The following object(s) are masked from 'package:base':

det

> a = Matrix((1:16)^2, nrow=4, ncol=4)
> a
4 x 4 Matrix of class "dgeMatrix"
[,1] [,2] [,3] [,4]
[1,]    1   25   81  169
[2,]    4   36  100  196
[3,]    9   49  121  225
[4,]   16   64  144  256
> lu = lu(a)
> lu
'MatrixFactorization' of Formal class 'denseLU' [package "Matrix"] with 3 slots
..@ x   : num [1:16] 16 0.0625 0.5625 0.25 64 ...
..@ perm: int [1:4] 4 4 3 4
..@ Dim : int [1:2] 4 4
> lue = expand(lu)
> lue
$L 4 x 4 Matrix of class "dtrMatrix" (unitriangular) [,1] [,2] [,3] [,4] [1,] 1.0000000 . . . [2,] 0.0625000 1.0000000 . . [3,] 0.5625000 0.6190476 1.0000000 . [4,] 0.2500000 0.9523810 1.0000000 1.0000000$U
4 x 4 Matrix of class "dtrMatrix"
[,1]          [,2]          [,3]          [,4]
[1,]  1.600000e+01  6.400000e+01  1.440000e+02  2.560000e+02
[2,]             .  2.100000e+01  7.200000e+01  1.530000e+02
[3,]             .             . -4.571429e+00 -1.371429e+01
[4,]             .             .             . -6.090614e-15

$P 4 x 4 sparse Matrix of class "pMatrix" [1,] . | . . [2,] . . . | [3,] . . | . [4,] | . . . > PLU = lue$P %*% lue$L %*% lue$U
> PLU
4 x 4 Matrix of class "dgeMatrix"
[,1] [,2] [,3] [,4]
[1,]    1   25   81  169
[2,]    4   36  100  196
[3,]    9   49  121  225
[4,]   16   64  144  256


# 範例二

> library(Matrix)
> x=matrix(rnorm(9),c(3,3)); x
[,1] [,2] [,3]
[1,] -0.6334882 -0.3915563 0.4906192
[2,] 0.4591368 0.5246114 0.6949097
[3,] -0.4435543 -1.5035618 -0.0191876
# 根据例子, 需要将矩阵转换为CsparseMatrix 类. why???
> lu(x)

unable to find an inherited method for function "lu", for signature "matrix"
> A = as(x,"CsparseMatrix")
> p=lu(A)
# 结果是'MatrixFactorization' of Formal class 'sparseLU'
> p
'MatrixFactorization' of Formal class 'sparseLU' [package "Matrix"] with 5 slots
..@ L :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
.. .. ..@ i : int [1:6] 0 1 2 1 2 2
.. .. ..@ p : int [1:4] 0 3 5 6
.. .. ..@ Dim : int [1:2] 3 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$: NULL .. .. .. ..$ : NULL
.. .. ..@ x : num [1:6] 1.000 0.700 -0.725 1.000 -0.196 ...
.. .. ..@ uplo : chr "L"
.. .. ..@ diag : chr "N"
..@ U :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
.. .. ..@ i : int [1:6] 0 0 1 0 1 2
.. .. ..@ p : int [1:4] 0 1 3 6
.. .. ..@ Dim : int [1:2] 3 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$: NULL .. .. .. ..$ : NULL
.. .. ..@ x : num [1:6] -0.633 -0.392 -1.229 0.491 -0.363 ...
.. .. ..@ uplo : chr "U"
.. .. ..@ diag : chr "N"
..@ p : int [1:3] 0 2 1
..@ q : int [1:3] 0 1 2
..@ Dim: int(0)
> p@L
3 x 3 sparse Matrix of class "dtCMatrix"
[1,] 1.0000000 . .
[2,] 0.7001776 1.0000000 .
[3,] -0.7247755 -0.1958845 1
> p@U
3 x 3 sparse Matrix of class "dtCMatrix"
[1,] -0.6334882 -0.3915563 0.4906192
[2,] . -1.2294029 -0.3627082
[3,] . . 0.9794496
# L*U 既得原来的矩阵, 行顺序可能不同
> p@L %*% p@U
3 x 3 sparse Matrix of class "dgCMatrix"
[1,] -0.6334882 -0.3915563 0.4906192
[2,] -0.4435543 -1.5035618 -0.0191876
[3,] 0.4591368 0.5246114 0.6949097
> A
3 x 3 sparse Matrix of class "dgCMatrix"
[1,] -0.6334882 -0.3915563 0.4906192
[2,] 0.4591368 0.5246114 0.6949097
[3,] -0.4435543 -1.5035618 -0.0191876


page revision: 3, last edited: 26 Oct 2011 04:50