elementary linear algebra in R-language

Educated neither in engineering nor economics, linear algebra is not familiar to me. Before opening matrix textbook, I try to preview what to do. Let's start from 2-dimensional examples.

Crane-Tortoise Arithmetic

Linear equations with two unknowns is known as "Crane-Tortoise Arithmetic" in Japanese elementary schools. It may sound odd, suppose total number of heads are three and feet are ten for a group consisting unknown number of cranes and tortoises.
R-function c() generates vector, and rbind() row-wise bounds vectors into a matrix.

> b = c(3,10) ; b
[1]  3 10
> head = c(1,1)
> feet = c(2,4)
> m = rbind(head,feet); m
     [,1] [,2]
head    1    1
feet    2    4
> solve(m, b)
[1] 1 2

The solver answers that the group consists of a crane and two tortoises. You can see "solve.R" in R-language source tarball, solve() usually calls LAPACK's LA_GSEV (gaussian elimination using LU decomposition), and R discards pivoting results.

Iterative methods: Jacobi

decompose the coefficient matrix m into diagonal one d and remainder r, inverse d, and iterate. production of matrices is %*% operator in R-language.

> b = c(3,10) ; b
[1]  3 10
> d_1 = c(1,0)
> d_2 = c(0,4)
> r_1 = c(0,1)
> r_2 = c(2,0)
> d = rbind(d_1,d_2); d
d
    [,1] [,2]
d_1    1    0
d_2    0    4
> r = rbind(r_1,r_2); r
    [,1] [,2]
r_1    0    1
r_2    2    0
> d_inv = solve(d); d_inv
     d_1  d_2
[1,]   1 0.00
[2,]   0 0.25
> i = c(0,0)
> x_ = i
> x = d_inv %*% (b - r %*% x_) ; x_ = x ; x
     [,1]
[1,]  3.0
[2,]  2.5
> x = d_inv %*% (b - r %*% x_) ; x_ = x ; x
     [,1]
[1,]  0.5
[2,]  1.0

[1,] 2.00
[2,] 2.25

[1,] 0.75
[2,] 1.50

[1,] 1.500
[2,] 2.125

[1,] 0.875
[2,] 1.750

[1,] 1.2500
[2,] 2.0625

[1,] 0.9375
[2,] 1.8750

[1,] 1.12500
[2,] 2.03125

[1,] 0.96875
[2,] 1.93750

[1,] 1.062500
[2,] 2.015625

[1,] 0.984375
[2,] 1.968750

[1,] 1.031250
[2,] 2.007812

[1,] 0.9921875
[2,] 1.9843750

[1,] 1.015625
[2,] 2.003906

[1,] 0.9960938
[2,] 1.9921875

[1,] 1.007812
[2,] 2.001953

[1,] 0.9980469
[2,] 1.9960938

[1,] 1.003906
[2,] 2.000977

[1,] 0.9990234
[2,] 1.9980469

iteration started from initial vector (0,0) eventually converses into (1,2).

More iterative method: Gauss-Seidel

Similar to Jacobi method, but in every stages perform product-sum operations row-wise, and computation of next row uses the result of previous row. As fixpoint of mappings:

• Jacobi: x = f(x)
• Gauss-Seidel: x = f_1(f_2(..(f_n(x)..))

I don't know why converge differently.
note that previous value x_ is needless and = should be written as <- in R-language standard.

> b = c(3,10) ; b
[1]  3 10
> x_ = 0 ; y_ = 0
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 3 1
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 2.0 1.5
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.50 1.75
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.250 1.875
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.1250 1.9375
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.06250 1.96875
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.031250 1.984375
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.015625 1.992188
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.007812 1.996094
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.003906 1.998047
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.001953 1.999023
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.000977 1.999512
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.000488 1.999756
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.000244 1.999878

Cramer's rule, or use of determinants

substitute columns in the coefficient matrix m with vector b and obtain determinants for resulting matrices. Finally divide with det(m).

> b = c(3,10) ; b
[1]  3 10
> m = rbind(c(1,1),c(2,4)) ; m
     [,1] [,2]
[1,]    1    1
[2,]    2    4
> det(m)
[1] 2
> m_1 = cbind(b,c(1,4)) ; m_1
      b 
[1,]  3 1
[2,] 10 4
> det(m_1)
[1] 2
> m_2 = cbind(c(1,2), b) ; m_2
        b
[1,] 1  3
[2,] 2 10
> det(m_2)
[1] 4
> c(det(m_1) / det(m), det(m_2) / det(m))
[1] 1 2

Goodies for myself

• listing defined objects with ls()
• removing defined object with rm(object)
• eigenvalues and so on

> m
     [,1] [,2]
[1,]    1    1
[2,]    2    4
> eigen(m)
$values
[1] 4.5615528 0.4384472

$vectors
           [,1]       [,2]
[1,] -0.2703230 -0.8719282
[2,] -0.9627697  0.4896337

> a = eigen(m)$vectors[,1]
> m %*% a
          [,1]
[1,] -1.233093
[2,] -4.391725
> eigen(m)$values[1] * eigen(m)$vectors[,1]
[1] -1.233093 -4.391725
> dim(m)
[1] 2 2
> rank(m)
[1] 1.5 3.0 1.5 4.0
> rownames(m) = c("upper", "lower")
> m
     [,1] [,2]
upper    1    1
lower    2    4
> colnames(m) = c("right", "left")
> m
     right left
upper    1      1
lower    2      4
> "this"
[1] "this"
> c("this")
[1] "this"
> m
     right left
upper    1      1
lower    2      4
> m[1,1] = 5
> m
     right left
upper    5      1
lower    2      4
> m[1,1]
[1] 5
> m[2,2]
[1] 4
> m[1]
[1] 5
> m[4]
[1] 4
> length(m)
[1] 4
> t <- matrix(c(1,5,3,7),ncol=2,byrow=TRUE)
> t
     [,1] [,2]
[1,]    1    5
[2,]    3    7
> class(t)
[1] "matrix"
> t <- as.table(t)
> class(t)
[1] "table"
> t
    A   B
A   1   5
B   3   7
> attributes(t)
$dim
[1] 2 2

$dimnames
$dimnames[[1]]
[1] "A" "B"

$dimnames[[2]]
[1] "A" "B"


$class
[1] "table"