---
title: "11 - Linear algebra in `caracas`"
author: Mikkel Meyer Andersen and Søren Højsgaard
date: "`r date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{11 - Linear algebra in `caracas`}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Introduction

```{r, message=FALSE, echo=FALSE}
library(caracas)
##packageVersion("caracas")
```

```{r, include = FALSE}
inline_code <- function(x) {
  x
}

if (!has_sympy()) {
  # SymPy not available, so the chunks shall not be evaluated
  knitr::opts_chunk$set(eval = FALSE)
  
  inline_code <- function(x) {
    deparse(substitute(x))
  }
}
```

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].

# Elementary matrix operations

## Creating matrices / vectors

We now show different ways to create a symbolic matrix:

```{r}
A <- matrix(c("a", "b", "0", "1"), 2, 2) %>% as_sym()
A
A <- matrix_(c("a", "b", "0", "1"), 2, 2) # note the '_' postfix
A
A <- as_sym("[[a, 0], [b, 1]]")
A

A2 <- matrix_(c("a", "b", "c", "1"), 2, 2)
A2

B <- matrix_(c("a", "b", "0", 
              "c", "c", "a"), 2, 3)
B

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

D <- diag_(c("a", "b")) # note the '_' postfix
D
```

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


## Matrix-matrix sum and product

```{r}
A + A2
A %*% B
```

## Hadamard (elementwise) product

```{r}
A * A2
```

## Vector operations

```{r}
x <- as_sym(paste0("x", 1:3))
x
x + x
1 / x
x / x
```


## Reciprocal matrix

```{r}
reciprocal_matrix(A2)
reciprocal_matrix(A2, num = "2*a")
```


## Matrix inverse; solve system of linear equations

Solve $Ax=b$ for $x$:

```{r}
inv(A)
x <- solve_lin(A, b)
x
A %*% x ## Sanity check
```

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

```{r}
M <- as_sym("[[1, 2, 3], [4, 5, 6]]")
pinv(M)
B <- as_sym("[[7], [8]]") 
B
z <- do_la(M, "pinv_solve", B)
print(z, rowvec = FALSE) # Do not print column vectors as transposed row vectors
```


# 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.


## QR decomposition

```{r}
A <- matrix(c("a", "0", "0", "1"), 2, 2) %>% as_sym()
A
qr_res <- QRdecomposition(A)
qr_res$Q
qr_res$R
```

## Eigenvalues and eigenvectors

```{r}
eigenval(A)
```

```{r}
evec <- eigenvec(A)
evec
evec1 <- evec[[1]]$eigvec
evec1
simplify(evec1)

lapply(evec, function(l) simplify(l$eigvec))
```


## Inverse, Penrose-Moore pseudo inverse

```{r}
inv(A)
pinv(cbind(A, A)) # pseudo inverse
```


## Additional functionality for linear algebra

`do_la` short for "do linear algebra"


```{r}
args(do_la)
```

The above functions can be called:

```{r}
do_la(A, "QRdecomposition") # == QRdecomposition(A)
do_la(A, "inv")             # == inv(A)
do_la(A, "eigenvec")        # == eigenvec(A)
do_la(A, "eigenvals")       # == eigenval(A)
```

### Characteristic polynomial

```{r}
cp <- do_la(A, "charpoly")
cp
as_expr(cp)
```

### Rank

```{r}
do_la(A, "rank")
```

### Cofactor

```{r}
A <- matrix(c("a", "b", "0", "1"), 2, 2) %>% as_sym()
A
do_la(A, "cofactor", 0, 1)
do_la(A, "cofactor_matrix")
```

### Echelon form

```{r}
do_la(cbind(A, A), "echelon_form")
```

### Cholesky factorisation

```{r}
B <- as_sym("[[9, 3*I], [-3*I, 5]]")
B
do_la(B, "cholesky")
```

### Gram Schmidt

```{r}
B <- t(as_sym("[[ 2, 3, 5 ], [3, 6, 2], [8, 3, 6]]"))
do_la(B, "GramSchmidt")
```


### Reduced row-echelon form (rref)

```{r}
B <- t(as_sym("[[ 2, 3, 5 ], [4, 6, 10], [8, 3, 6] ]"))
B
B_rref <- do_la(B, "rref")
B_rref
```

### Column space, row space and null space

```{r}
B <- matrix(c(1, 3, 0, -2, -6, 0, 3, 9, 6), nrow = 3) %>% as_sym()
B
columnspace(B)
rowspace(B)
x <- nullspace(B)
x
rref(B)
B %*% x
```

### Singular values, svd

```{r}
B <- t(as_sym("[[ 2, 3, 5 ], [4, 6, 10], [8, 3, 6], [8, 3, 6] ]"))
B
singular_values(B)
```


<!-- ```{r} -->
<!-- ##svd(B) FIXME -->
<!-- # # FIXME -->
<!-- # B <- t(as_sym("[[ 2, 3, 5 ], [4, 6, 10], [8, 3, 6], [8, 3, 6] ]")) -->
<!-- # B -->
<!-- #  -->
<!-- # #do_la(B %*% t(B), "left_eigenvects") -->
<!-- #  -->
<!-- # U <- eigenvec(B %*% t(B)) %>% lapply(function(v) v$eigvec) %>% do.call(cbind, .) -->
<!-- # V <- eigenvec(t(B) %*% B) %>% lapply(function(v) v$eigvec) %>% do.call(cbind, .) -->
<!-- #  -->
<!-- # d <- singular_values(B)[seq_len(ncol(U))] -->
<!-- # D <- diag(ncol(U)) %>% as_sym() -->
<!-- # for (i in seq_along(d)) { -->
<!-- #   D[i, i] <- as.character(d[[i]]) -->
<!-- # } -->
<!-- #  -->
<!-- # d <- singular_values(B) %>% do.call(cbind, .) -->
<!-- # D <- as_diag(d) -->
<!-- #  -->
<!-- #  -->
<!-- # r_svd <- as_expr(B) %>% svd() -->
<!-- #  -->
<!-- # r_svd$d -->
<!-- # d %>% as_expr() -->
<!-- #  -->
<!-- # r_svd$u -->
<!-- # U %>% as_expr() -->
<!-- #  -->
<!-- # r_svd$v -->
<!-- # V %>% as_expr() -->
<!-- #  -->
<!-- # r_svd$u %*% diag(r_svd$d) %*% t(r_svd$v) -->
<!-- # U %*% D %*% t(V) %>% as_expr() -->
<!-- ``` -->