---
title: "Cramer's rule in civilised form with Clifford algebra"
author: "Robin K. S. Hankin"
bibliography: clifford.bib
vignette: >
  %\VignetteIndexEntry{Cramer's rule using Clifford algebra}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r out.width='20%', out.extra='style="float:right; padding:10px"',echo=FALSE}
knitr::include_graphics(system.file("help/figures/clifford.png", package = "clifford"))
```

To cite the `clifford` package in publications please use
@hankin2022_clifford.  This short document shows a nice application of
Clifford algebras to linear algebra.  Suppose we have vectors
${\mathbf a}, {\mathbf b}, {\mathbf c}$ spanning $\mathbb{R}^3$ and
are given ${\mathbf x}\in\mathbb{R}^3$.  We wish to write ${\mathbf
x}=\alpha {\mathbf a}+\beta {\mathbf b}+\gamma {\mathbf c}$ for some
$\alpha,\beta,\gamma\in\mathbb{R}$.  The traditional Cramer's rule for
finding $\alpha,\beta,\gamma$ would be

\[
\alpha=\frac{
	\det\begin{bmatrix}
	x_1&b_1&c_1\\
	x_2&b_2&c_2\\
	x_3&b_3&c_3
	\end{bmatrix}
	}{
	\det\begin{bmatrix}
	a_1&b_1&c_1\\
	a_2&b_2&c_2\\
	a_3&b_3&c_3
	\end{bmatrix}
}
\qquad\beta=\frac{
	\det\begin{bmatrix}
	a_1&x_1&c_1\\
	a_2&x_2&c_2\\
	a_3&x_3&c_3
	\end{bmatrix}
	}{
	\det\begin{bmatrix}
	a_1&b_1&c_1\\
	a_2&b_2&c_2\\
	a_3&b_3&c_3
	\end{bmatrix}
}
\qquad
\gamma=\frac{
	\det\begin{bmatrix}
	a_1&b_1&x_1\\
	a_2&b_2&x_2\\
	a_3&b_3&x_3
	\end{bmatrix}
	}{
	\det\begin{bmatrix}
	a_1&b_1&c_1\\
	a_2&b_2&c_2\\
	a_3&b_3&c_3
	\end{bmatrix}
}
\]

where ${\mathbf x}=(x_1,x_2,x_3)^T$, ${\mathbf a}=(a_1,a_2,a_3)^T$,
${\mathbf b}=(b_1,b_2,b_3)^T$ and ${\mathbf c}=(c_1,c_2,c_3)^T$.
However, observe that this solution, while accurate, requires one to
take a coordinate basis; and offers little in the way of intuition.

## Using Clifford algebra

Considering $\mathbb{R}^3$ as a vector space and given vectors
${\mathbf a}, {\mathbf b}, {\mathbf c}$ spanning the space we can
express any vector ${\mathbf x}\in\mathbb{R}^3$ as

\[\mathbf{x}=
\left(\frac{{\mathbf x}\wedge{\mathbf b}\wedge{\mathbf c}}{{\mathbf a}\wedge{\mathbf b}\wedge{\mathbf c}}\right){\mathbf a}+
\left(\frac{{\mathbf a}\wedge{\mathbf x}\wedge{\mathbf c}}{{\mathbf a}\wedge{\mathbf b}\wedge{\mathbf c}}\right){\mathbf b}+
\left(\frac{{\mathbf a}\wedge{\mathbf b}\wedge{\mathbf x}}{{\mathbf a}\wedge{\mathbf b}\wedge{\mathbf c}}\right){\mathbf c}
\]

which is Cramer's rule expressed directly in vector form (rather than
components).  Observe that the numerator and denominator of each
bracketed term is a pseudoscalar; the ratio of two pseudoscalars is an
ordinary scalar.  Package idiom is straightforward:

```{r, label=loadlib,results="hide",include=FALSE}
library("clifford",quietly=TRUE)  # document requires package version 1.0-9 or above
set.seed(0)
```

```{r label=useone,eval=TRUE}
a <- as.1vector(runif(3))
b <- as.1vector(runif(3))
c <- as.1vector(runif(3))

(x <- as.1vector(1:3))

options(maxdim = 3)  # needed to drop() pseudoscalars
abc <- drop(a ^ b ^ c)

alpha <- drop(x ^ b ^ c)/abc
beta  <- drop(a ^ x ^ c)/abc
gamma <- drop(a ^ b ^ x)/abc

c(alpha,beta,gamma)

alpha*a + beta*b + gamma*c
Mod(alpha*a + beta*b + gamma*c-x)
```

Thus we have expressed ${\mathbf x}$ (except for possible roundoff
error) as a linear combination of ${\mathbf a},{\mathbf b},{\mathbf
c}$, specifically ${\mathbf x}=\alpha{\mathbf a}+\beta{\mathbf
b}+\gamma{\mathbf c}$.  Conversely, we might know the coefficients and
try to determine them using package idiom.  Here we will use $1,2,3$
and suppose that ${\mathbf y}=1{\mathbf a}+2{\mathbf b}+3{\mathbf c}$:

```{r label=showdrop}
y <- a*1 + b*2 + c*3
c(
    drop(y ^ b ^ c)/abc,
    drop(a ^ y ^ c)/abc,
    drop(a ^ b ^ y)/abc
)
```


# Higher dimensional space

To accomplish this in arbitrary-dimensional space is straightforward.
Here we consider $\mathbb{R}^{5}$:

```{r arbdim}
n <- 5                                                 # dimensionality of space
options(maxdim=5)                                      # safety precaution
x <- as.1vector(seq_len(n))                            # target vector
x
L <- replicate(n,as.1vector(rnorm(n)),simplify=FALSE)  # spanning vectors
subst <- function(L,n,x){L[[n]] <- x; return(L)}       # list substitution
coeff <- function(n,L,x){
      drop(Reduce(`^`,subst(L,n,x))/Reduce(`^`,L))
}
```

Then the coefficients are given by:

```{r showcoeffs}
(alpha <- sapply(seq_len(n),coeff,L,x))
```

and we can reconstitute vector $x$:

```{r reconvec}
out <- as.clifford(0)
f <- function(i){alpha[i]*L[[i]]}
for(i in seq_len(n)){
      out <- out + f(i)
}
Mod(out-x)  # zero to numerical precision
```

Or, somewhat slicker:

```{r somewhatslicker}
Reduce(`+`,sapply(seq_len(n),f,simplify=FALSE))
```


Conversely, if we know the coefficients are, say, `15:11`, then we would have

```{r conversely}
coeffs <- 15:11
x <- 0
for(i in seq_len(5)){x <- x + coeffs[i]*L[[i]]}
x
```

And then to find the coefficients:

```{r findcoeffs}
sapply(seq_len(n),coeff,L,x)
```

Above we see that the original coefficients are recovered, up to numerical accuracy.

## References

```{r restoredefaults, echo=FALSE}
options(maxdim=NULL) # restore default for maxdim: options persist
                     # between vignettes, and leaving maxdim set
                     # causes problems for other vignettes
```