The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

11 - Linear algebra in caracas

Mikkel Meyer Andersen and Søren Højsgaard

Thu Nov 30 15:39:50 2023

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
#> c: ⎡a  0⎤
#>    ⎢    ⎥
#>    ⎣b  1⎦
A <- matrix_(c("a", "b", "0", "1"), 2, 2) # note the '_' postfix
A
#> c: ⎡a  0⎤
#>    ⎢    ⎥
#>    ⎣b  1⎦
A <- as_sym("[[a, 0], [b, 1]]")
A
#> c: ⎡a  0⎤
#>    ⎢    ⎥
#>    ⎣b  1⎦

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

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

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

D <- diag_(c("a", "b")) # note the '_' postfix
D
#> c: ⎡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
#> c: ⎡2⋅a  c⎤
#>    ⎢      ⎥
#>    ⎣2⋅b  2⎦
A %*% B
#> c: ⎡   2               ⎤
#>    ⎢  a      0    a⋅c  ⎥
#>    ⎢                   ⎥
#>    ⎣a⋅b + b  c  a + b⋅c⎦

2.3 Hadamard (elementwise) product

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

2.4 Vector operations

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

2.5 Reciprocal matrix

reciprocal_matrix(A2)
#> c: ⎡1  1⎤
#>    ⎢─  ─⎥
#>    ⎢a  c⎥
#>    ⎢    ⎥
#>    ⎢1   ⎥
#>    ⎢─  1⎥
#>    ⎣b   ⎦
reciprocal_matrix(A2, num = "2*a")
#> c: ⎡     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)
#> c: ⎡ 1    ⎤
#>    ⎢ ─   0⎥
#>    ⎢ a    ⎥
#>    ⎢      ⎥
#>    ⎢-b    ⎥
#>    ⎢───  1⎥
#>    ⎣ a    ⎦
x <- solve_lin(A, b)
x
#> c: ⎡b₁       b⋅b₁⎤
#>    ⎢──  b₂ - ────⎥
#>    ⎣a         a  ⎦ᵀ
A %*% x ## Sanity check
#> c: [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)
#> c: ⎡-17       ⎤
#>    ⎢────  4/9 ⎥
#>    ⎢ 18       ⎥
#>    ⎢          ⎥
#>    ⎢-1/9  1/9 ⎥
#>    ⎢          ⎥
#>    ⎢ 13       ⎥
#>    ⎢ ──   -2/9⎥
#>    ⎣ 18       ⎦
B <- as_sym("[[7], [8]]") 
B
#> c: [7  8]ᵀ
z <- do_la(M, "pinv_solve", B)
print(z, rowvec = FALSE) # Do not print column vectors as transposed row vectors
#> c: ⎡ 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
#> c: ⎡a  0⎤
#>    ⎢    ⎥
#>    ⎣0  1⎦
qr_res <- QRdecomposition(A)
qr_res$Q
#> c: ⎡ a    ⎤
#>    ⎢───  0⎥
#>    ⎢│a│   ⎥
#>    ⎢      ⎥
#>    ⎣ 0   1⎦
qr_res$R
#> c: ⎡│a│  0⎤
#>    ⎢      ⎥
#>    ⎣ 0   1⎦

3.2 Eigenvalues and eigenvectors

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

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

3.3 Inverse, Penrose-Moore pseudo inverse

inv(A)
#> c: ⎡1   ⎤
#>    ⎢─  0⎥
#>    ⎢a   ⎥
#>    ⎢    ⎥
#>    ⎣0  1⎦
pinv(cbind(A, A)) # pseudo inverse
#> c: ⎡ 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
#> c: ⎡ a    ⎤
#>    ⎢───  0⎥
#>    ⎢│a│   ⎥
#>    ⎢      ⎥
#>    ⎣ 0   1⎦
#> 
#> $R
#> c: ⎡│a│  0⎤
#>    ⎢      ⎥
#>    ⎣ 0   1⎦
do_la(A, "inv")             # == inv(A)
#> c: ⎡1   ⎤
#>    ⎢─  0⎥
#>    ⎢a   ⎥
#>    ⎢    ⎥
#>    ⎣0  1⎦
do_la(A, "eigenvec")        # == eigenvec(A)
#> [[1]]
#> [[1]]$eigval
#> c: 1
#> 
#> [[1]]$eigmult
#> [1] 1
#> 
#> [[1]]$eigvec
#> c: [0  1]ᵀ
#> 
#> 
#> [[2]]
#> [[2]]$eigval
#> c: a
#> 
#> [[2]]$eigmult
#> [1] 1
#> 
#> [[2]]$eigvec
#> c: [1  0]ᵀ
do_la(A, "eigenvals")       # == eigenval(A)
#> [[1]]
#> [[1]]$eigval
#> c: a
#> 
#> [[1]]$eigmult
#> [1] 1
#> 
#> 
#> [[2]]
#> [[2]]$eigval
#> c: 1
#> 
#> [[2]]$eigmult
#> [1] 1

3.4.1 Characteristic polynomial

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

3.4.2 Rank

do_la(A, "rank")
#> c: 2

3.4.3 Cofactor

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

3.4.4 Echelon form

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

3.4.5 Cholesky factorisation

B <- as_sym("[[9, 3*I], [-3*I, 5]]")
B
#> c: ⎡ 9    3⋅ⅈ⎤
#>    ⎢         ⎥
#>    ⎣-3⋅ⅈ   5 ⎦
do_la(B, "cholesky")
#> c: ⎡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")
#> c: ⎡    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
#> c: ⎡2  4   8⎤
#>    ⎢        ⎥
#>    ⎢3  6   3⎥
#>    ⎢        ⎥
#>    ⎣5  10  6⎦
B_rref <- do_la(B, "rref")
B_rref
#> $mat
#> c: ⎡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
#> c: ⎡1  -2  3⎤
#>    ⎢        ⎥
#>    ⎢3  -6  9⎥
#>    ⎢        ⎥
#>    ⎣0  0   6⎦
columnspace(B)
#> c: ⎡1  3⎤
#>    ⎢    ⎥
#>    ⎢3  9⎥
#>    ⎢    ⎥
#>    ⎣0  6⎦
rowspace(B)
#> c: [1  -2  3  0  0  6]
x <- nullspace(B)
x
#> c: [2  1  0]ᵀ
rref(B)
#> $mat
#> c: ⎡1  -2  0⎤
#>    ⎢        ⎥
#>    ⎢0  0   1⎥
#>    ⎢        ⎥
#>    ⎣0  0   0⎦
#> 
#> $pivot_vars
#> [1] 1 3
B %*% x
#> c: [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
#> c: ⎡2  4   8  8⎤
#>    ⎢           ⎥
#>    ⎢3  6   3  3⎥
#>    ⎢           ⎥
#>    ⎣5  10  6  6⎦
singular_values(B)
#> [[1]]
#> c:   ______________
#>    ╲╱ √30446 + 204
#> 
#> [[2]]
#> c:   ______________
#>    ╲╱ 204 - √30446
#> 
#> [[3]]
#> c: 0
#> 
#> [[4]]
#> c: 0

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.