Linear algebra in caracas

Mikkel Meyer Andersen and Søren Højsgaard

Fri Feb 11 08:32:16 2022

1 Introduction

This vignette is based on caracas version r packageVersion("caracas"). caracas is avavailable on CRAN at [https://cran.r-project.org/package=caracas] and on github at [https://github.com/r-cas/caracas].

2 Elementary matrix operations

2.1 Creating matrices / vectors

We now show different ways to create a symbolic matrix:

A <- matrix(c("a", "b", "0", "1"), 2, 2) %>% as_sym()
A
#> [caracas]: ⎡a  0⎤
#>            ⎢    ⎥
#>            ⎣b  1⎦
A <- matrix_(c("a", "b", "0", "1"), 2, 2) # note the '_' postfix
A
#> [caracas]: ⎡a  0⎤
#>            ⎢    ⎥
#>            ⎣b  1⎦
A <- as_sym("[[a, 0], [b, 1]]")
A
#> [caracas]: ⎡a  0⎤
#>            ⎢    ⎥
#>            ⎣b  1⎦

A2 <- matrix_(c("a", "b", "c", "1"), 2, 2)
A2
#> [caracas]: ⎡a  c⎤
#>            ⎢    ⎥
#>            ⎣b  1⎦

B <- matrix_(c("a", "b", "0", 
              "c", "c", "a"), 2, 3)
B
#> [caracas]: ⎡a  0  c⎤
#>            ⎢       ⎥
#>            ⎣b  c  a⎦

b <- matrix_(c("b1", "b2"), nrow = 2)

D <- diag_(c("a", "b")) # note the '_' postfix
D
#> [caracas]: ⎡a  0⎤
#>            ⎢    ⎥
#>            ⎣0  b⎦

Note that a vector is a matrix in which one of the dimensions is one.

2.2 Matrix-matrix sum and product

A + A2
#> [caracas]: ⎡2⋅a  c⎤
#>            ⎢      ⎥
#>            ⎣2⋅b  2⎦
A %*% B
#> [caracas]: ⎡   2               ⎤
#>            ⎢  a      0    a⋅c  ⎥
#>            ⎢                   ⎥
#>            ⎣a⋅b + b  c  a + b⋅c⎦

2.3 Hadamard (elementwise) product

A * A2
#> [caracas]: ⎡ 2   ⎤
#>            ⎢a   0⎥
#>            ⎢     ⎥
#>            ⎢ 2   ⎥
#>            ⎣b   1⎦

2.4 Vector operations

x <- as_sym(paste0("x", 1:3))
x
#> [caracas]: [x₁  x₂  x₃]ᵀ
x + x
#> [caracas]: [2⋅x₁  2⋅x₂  2⋅x₃]ᵀ
1 / x
#> [caracas]: ⎡1   1   1 ⎤
#>            ⎢──  ──  ──⎥
#>            ⎣x₁  x₂  x₃⎦ᵀ
x / x
#> [caracas]: [1  1  1]ᵀ

2.5 Reciprocal matrix

reciprocal_matrix(A2)
#> [caracas]: ⎡1  1⎤
#>            ⎢─  ─⎥
#>            ⎢a  c⎥
#>            ⎢    ⎥
#>            ⎢1   ⎥
#>            ⎢─  1⎥
#>            ⎣b   ⎦
reciprocal_matrix(A2, num = "2*a")
#> [caracas]: ⎡     2⋅a⎤
#>            ⎢ 2   ───⎥
#>            ⎢      c ⎥
#>            ⎢        ⎥
#>            ⎢2⋅a     ⎥
#>            ⎢───  2⋅a⎥
#>            ⎣ b      ⎦

2.6 Matrix inverse; solve system of linear equations

Solve \(Ax=b\) for \(x\):

inv(A)
#> [caracas]: ⎡ 1    ⎤
#>            ⎢ ─   0⎥
#>            ⎢ a    ⎥
#>            ⎢      ⎥
#>            ⎢-b    ⎥
#>            ⎢───  1⎥
#>            ⎣ a    ⎦
x <- solve_lin(A, b)
x
#> [caracas]: ⎡b₁       b⋅b₁⎤
#>            ⎢──  b₂ - ────⎥
#>            ⎣a         a  ⎦ᵀ
A %*% x ## Sanity check
#> [caracas]: [b₁  b₂]ᵀ

2.7 Generalized (Penrose-Moore) inverse; solve system of linear equations [TBW]

M <- as_sym("[[1, 2, 3], [4, 5, 6]]")
pinv(M)
#> [caracas]: ⎡-17       ⎤
#>            ⎢────  4/9 ⎥
#>            ⎢ 18       ⎥
#>            ⎢          ⎥
#>            ⎢-1/9  1/9 ⎥
#>            ⎢          ⎥
#>            ⎢ 13       ⎥
#>            ⎢ ──   -2/9⎥
#>            ⎣ 18       ⎦
B <- as_sym("[[7], [8]]") 
B
#> [caracas]: [7  8]ᵀ
z <- do_la(M, "pinv_solve", B)
print(z, rowvec = FALSE) # Do not print column vectors as transposed row vectors
#> [caracas]: ⎡ w₀ ₀   w₁ ₀   w₂ ₀   55  ⎤
#>            ⎢ ──── - ──── + ──── - ──  ⎥
#>            ⎢  6      3      6     18  ⎥
#>            ⎢                          ⎥
#>            ⎢  w₀ ₀   2⋅w₁ ₀   w₂ ₀   1⎥
#>            ⎢- ──── + ────── - ──── + ─⎥
#>            ⎢   3       3       3     9⎥
#>            ⎢                          ⎥
#>            ⎢ w₀ ₀   w₁ ₀   w₂ ₀   59  ⎥
#>            ⎢ ──── - ──── + ──── + ──  ⎥
#>            ⎣  6      3      6     18  ⎦

3 More special linear algebra functionality

Below we present convenient functions for performing linear algebra operations. If you need a function in SymPy for which we have not supplied a convinience function (see https://docs.sympy.org/latest/modules/matrices/matrices.html), you can still call it with the do_la() (short for “do linear algebra”) function presented at the end of this section.

3.1 QR decomposition

A <- matrix(c("a", "0", "0", "1"), 2, 2) %>% as_sym()
A
#> [caracas]: ⎡a  0⎤
#>            ⎢    ⎥
#>            ⎣0  1⎦
qr_res <- QRdecomposition(A)
qr_res$Q
#> [caracas]: ⎡ a    ⎤
#>            ⎢───  0⎥
#>            ⎢│a│   ⎥
#>            ⎢      ⎥
#>            ⎣ 0   1⎦
qr_res$R
#> [caracas]: ⎡│a│  0⎤
#>            ⎢      ⎥
#>            ⎣ 0   1⎦

3.2 Eigenvalues and eigenvectors

eigenval(A)
#> [[1]]
#> [[1]]$eigval
#> [caracas]: a
#> 
#> [[1]]$eigmult
#> [1] 1
#> 
#> 
#> [[2]]
#> [[2]]$eigval
#> [caracas]: 1
#> 
#> [[2]]$eigmult
#> [1] 1
evec <- eigenvec(A)
evec
#> [[1]]
#> [[1]]$eigval
#> [caracas]: 1
#> 
#> [[1]]$eigmult
#> [1] 1
#> 
#> [[1]]$eigvec
#> [caracas]: [0  1]ᵀ
#> 
#> 
#> [[2]]
#> [[2]]$eigval
#> [caracas]: a
#> 
#> [[2]]$eigmult
#> [1] 1
#> 
#> [[2]]$eigvec
#> [caracas]: [1  0]ᵀ
evec1 <- evec[[1]]$eigvec
evec1
#> [caracas]: [0  1]ᵀ
simplify(evec1)
#> [caracas]: [0  1]ᵀ

lapply(evec, function(l) simplify(l$eigvec))
#> [[1]]
#> [caracas]: [0  1]ᵀ
#> 
#> [[2]]
#> [caracas]: [1  0]ᵀ

3.3 Inverse, Penrose-Moore pseudo inverse

inv(A)
#> [caracas]: ⎡1   ⎤
#>            ⎢─  0⎥
#>            ⎢a   ⎥
#>            ⎢    ⎥
#>            ⎣0  1⎦
pinv(cbind(A, A)) # pseudo inverse
#> [caracas]: ⎡ 1      ⎤
#>            ⎢───   0 ⎥
#>            ⎢2⋅a     ⎥
#>            ⎢        ⎥
#>            ⎢ 0   1/2⎥
#>            ⎢        ⎥
#>            ⎢ 1      ⎥
#>            ⎢───   0 ⎥
#>            ⎢2⋅a     ⎥
#>            ⎢        ⎥
#>            ⎣ 0   1/2⎦

3.4 Additional functionality for linear algebra

do_la short for “do linear algebra”

args(do_la)
#> function (x, slot, ...) 
#> NULL

The above functions can be called:

do_la(A, "QRdecomposition") # == QRdecomposition(A)
#> $Q
#> [caracas]: ⎡ a    ⎤
#>            ⎢───  0⎥
#>            ⎢│a│   ⎥
#>            ⎢      ⎥
#>            ⎣ 0   1⎦
#> 
#> $R
#> [caracas]: ⎡│a│  0⎤
#>            ⎢      ⎥
#>            ⎣ 0   1⎦
do_la(A, "inv")             # == inv(A)
#> [caracas]: ⎡1   ⎤
#>            ⎢─  0⎥
#>            ⎢a   ⎥
#>            ⎢    ⎥
#>            ⎣0  1⎦
do_la(A, "eigenvec")        # == eigenvec(A)
#> [[1]]
#> [[1]]$eigval
#> [caracas]: 1
#> 
#> [[1]]$eigmult
#> [1] 1
#> 
#> [[1]]$eigvec
#> [caracas]: [0  1]ᵀ
#> 
#> 
#> [[2]]
#> [[2]]$eigval
#> [caracas]: a
#> 
#> [[2]]$eigmult
#> [1] 1
#> 
#> [[2]]$eigvec
#> [caracas]: [1  0]ᵀ
do_la(A, "eigenvals")       # == eigenval(A)
#> [[1]]
#> [[1]]$eigval
#> [caracas]: a
#> 
#> [[1]]$eigmult
#> [1] 1
#> 
#> 
#> [[2]]
#> [[2]]$eigval
#> [caracas]: 1
#> 
#> [[2]]$eigmult
#> [1] 1

3.4.1 Characteristic polynomium

cp <- do_la(A, "charpoly")
cp
#> [caracas]:      2             
#>            a + λ  + λ⋅(-a - 1)
as_expr(cp)
#> expression(a + lambda^2 + lambda * (-a - 1))

3.4.2 Rank

do_la(A, "rank")
#> [caracas]: 2

3.4.3 Cofactor

A <- matrix(c("a", "b", "0", "1"), 2, 2) %>% as_sym()
A
#> [caracas]: ⎡a  0⎤
#>            ⎢    ⎥
#>            ⎣b  1⎦
do_la(A, "cofactor", 0, 1)
#> [caracas]: -1.0⋅b
do_la(A, "cofactor_matrix")
#> [caracas]: ⎡1  -b⎤
#>            ⎢     ⎥
#>            ⎣0  a ⎦

3.4.4 Echelon form

do_la(cbind(A, A), "echelon_form")
#> [caracas]: ⎡a  0  a  0⎤
#>            ⎢          ⎥
#>            ⎣0  a  0  a⎦

3.4.5 Cholesky factorisation

B <- as_sym("[[9, 3*I], [-3*I, 5]]")
B
#> [caracas]: ⎡ 9    3⋅ⅈ⎤
#>            ⎢         ⎥
#>            ⎣-3⋅ⅈ   5 ⎦
do_la(B, "cholesky")
#> [caracas]: ⎡3   0⎤
#>            ⎢     ⎥
#>            ⎣-ⅈ  2⎦

3.4.6 Gram Schmidt

B <- t(as_sym("[[ 2, 3, 5 ], [3, 6, 2], [8, 3, 6]]"))
do_la(B, "GramSchmidt")
#> [caracas]: ⎡    23    1692 ⎤
#>            ⎢2   ──    ──── ⎥
#>            ⎢    19    353  ⎥
#>            ⎢               ⎥
#>            ⎢    63   -1551 ⎥
#>            ⎢3   ──   ──────⎥
#>            ⎢    19    706  ⎥
#>            ⎢               ⎥
#>            ⎢   -47   -423  ⎥
#>            ⎢5  ────  ───── ⎥
#>            ⎣    19    706  ⎦

3.4.7 Reduced row-echelon form (rref)

B <- t(as_sym("[[ 2, 3, 5 ], [4, 6, 10], [8, 3, 6] ]"))
B
#> [caracas]: ⎡2  4   8⎤
#>            ⎢        ⎥
#>            ⎢3  6   3⎥
#>            ⎢        ⎥
#>            ⎣5  10  6⎦
B_rref <- do_la(B, "rref")
B_rref
#> $mat
#> [caracas]: ⎡1  2  0⎤
#>            ⎢       ⎥
#>            ⎢0  0  1⎥
#>            ⎢       ⎥
#>            ⎣0  0  0⎦
#> 
#> $pivot_vars
#> [1] 1 3

3.4.8 Column space, row space and null space

B <- matrix(c(1, 3, 0, -2, -6, 0, 3, 9, 6), nrow = 3) %>% as_sym()
B
#> [caracas]: ⎡1  -2  3⎤
#>            ⎢        ⎥
#>            ⎢3  -6  9⎥
#>            ⎢        ⎥
#>            ⎣0  0   6⎦
columnspace(B)
#> [[1]]
#> [caracas]: [1  3  0]ᵀ
#> 
#> [[2]]
#> [caracas]: [3  9  6]ᵀ
rowspace(B)
#> [[1]]
#> [caracas]: [1  -2  3]
#> 
#> [[2]]
#> [caracas]: [0  0  6]
x <- nullspace(B)
x
#> [[1]]
#> [caracas]: [2  1  0]ᵀ
rref(B)
#> $mat
#> [caracas]: ⎡1  -2  0⎤
#>            ⎢        ⎥
#>            ⎢0  0   1⎥
#>            ⎢        ⎥
#>            ⎣0  0   0⎦
#> 
#> $pivot_vars
#> [1] 1 3
B %*% x[[1]]
#> [caracas]: [0  0  0]ᵀ

3.4.9 Singular values, svd

B <- t(as_sym("[[ 2, 3, 5 ], [4, 6, 10], [8, 3, 6], [8, 3, 6] ]"))
B
#> [caracas]: ⎡2  4   8  8⎤
#>            ⎢           ⎥
#>            ⎢3  6   3  3⎥
#>            ⎢           ⎥
#>            ⎣5  10  6  6⎦
singular_values(B)
#> [[1]]
#> [caracas]:   ______________
#>            ╲╱ √30446 + 204
#> 
#> [[2]]
#> [caracas]:   ______________
#>            ╲╱ 204 - √30446
#> 
#> [[3]]
#> [caracas]: 0
#> 
#> [[4]]
#> [caracas]: 0