caracasThis 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].
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.
A + A2
#> [caracas]: ⎡2⋅a c⎤
#> ⎢ ⎥
#> ⎣2⋅b 2⎦
A %*% B
#> [caracas]: ⎡ 2 ⎤
#> ⎢ a 0 a⋅c ⎥
#> ⎢ ⎥
#> ⎣a⋅b + b c a + b⋅c⎦A * A2
#> [caracas]: ⎡ 2 ⎤
#> ⎢a 0⎥
#> ⎢ ⎥
#> ⎢ 2 ⎥
#> ⎣b 1⎦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]ᵀ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 ⎦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₂]ᵀ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 ⎦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.
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⎦eigenval(A)
#> [[1]]
#> [[1]]$eigval
#> [caracas]: a
#>
#> [[1]]$eigmult
#> [1] 1
#>
#>
#> [[2]]
#> [[2]]$eigval
#> [caracas]: 1
#>
#> [[2]]$eigmult
#> [1] 1evec <- 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]ᵀ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⎦do_la short for “do linear algebra”
args(do_la)
#> function (x, slot, ...)
#> NULLThe 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] 1cp <- do_la(A, "charpoly")
cp
#> [caracas]: 2
#> a + λ + λ⋅(-a - 1)
as_expr(cp)
#> expression(a + lambda^2 + lambda * (-a - 1))do_la(A, "rank")
#> [caracas]: 2A <- 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 ⎦do_la(cbind(A, A), "echelon_form")
#> [caracas]: ⎡a 0 a 0⎤
#> ⎢ ⎥
#> ⎣0 a 0 a⎦B <- as_sym("[[9, 3*I], [-3*I, 5]]")
B
#> [caracas]: ⎡ 9 3⋅ⅈ⎤
#> ⎢ ⎥
#> ⎣-3⋅ⅈ 5 ⎦
do_la(B, "cholesky")
#> [caracas]: ⎡3 0⎤
#> ⎢ ⎥
#> ⎣-ⅈ 2⎦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 ⎦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 3B <- 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]ᵀ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