MoA Interpretation of the Iterative Conjugate Gradient Method with Psi Reduction - A Tutorial to teach the Mathematically literate in Linear and Tensor Algebra: Part I

TL;DR

This tutorial explains the MoA interpretation of the iterative conjugate gradient method with Psi reduction, aiming to aid mathematically literate individuals in understanding advanced linear and tensor algebra concepts for applications in AI, HPC, and more.

Contribution

It provides a clear, evolving tutorial on MoA and Psi calculus, connecting foundational mathematics to practical computational methods and design verification.

Findings

Enhanced understanding of MoA interpretation for conjugate gradient methods

Clarification of Psi calculus applications in tensor algebra and algorithm design

Facilitation of learning for mathematically literate audiences in advanced linear algebra

Abstract

It is often difficult to learn new mathematics semantically and syntactically, even when there are similarities in the words and meaning when discussed aloud. The goal of this document is to facilitate learning through explanations and definitions relating our common mathematical knowledge and highlighting what is new. It is meant to be a working document that will evolve based on feedback from target audiences, those mathematically literate in linear and tensor algebra, those that want to learn MoA, Psi Calculus, and its uses, those that want and need the ability to prove a design, either in hardware or software through the ONF, Operational Normal Form, and those wanting to exploit all resources optimally, especially when Tensor Algebra, i.e. algorithms foundational to their application,are needed: Knowledge Representation, Machine Learning, Signal Processing, AI, HPC, etc.

Equations

$$A\vec{x}=\vec{b}$$

$$A_{+}{\bullet }_{\times }\vec{x}\equiv \vec{b}$$

$$\rho A\equiv <nn>$$

$$\rho \vec{x}\equiv \rho \vec{b}\equiv <n>$$

$${\mathbf{r}}_0$$

$${\mathbf{p}}_0$$

$$<0>\psi X$$

$$<0>\psi P$$

$$k:=0$$

$${\alpha }_k:=\frac{{\mathbf{r}}_k^T{\mathbf{r}}_k}{{\mathbf{p}}_k^T\mathbf{A}{\mathbf{p}}_k}$$

$${\mathbf{x}}_{k+1}:={\mathbf{x}}_k+{\alpha }_k{\mathbf{p}}_k$$

$${\mathbf{r}}_{k+1}:={\mathbf{r}}_k-{\alpha }_k\mathbf{A}{\mathbf{p}}_k$$

$${\mathbf{r}}_{k+1}$$

$${\beta }_k:=\frac{{\mathbf{r}}_{k+1}^T{\mathbf{r}}_{k+1}}{{\mathbf{r}}_k^T{\mathbf{r}}_k}$$

$${\mathbf{p}}_{k+1}:={\mathbf{r}}_{k+1}+{\beta }_k{\mathbf{p}}_k$$

$$k:=k+1$$

$${\mathbf{x}}_{k+1}$$

$$0\le i<n$$

$$<i>\psi \vec{\alpha }$$

$$<i+1>\psi X$$

$$<i+1>\psi R$$

$$<i>\psi \vec{\beta }$$

$$<i+1>\psi P$$

$$<i+1>\psi X$$

$$<i+1>\psi R$$

$$<i+1>\psi P$$

$$<i+1>\psi X$$

$$<i-i+1>\psi X$$

$$<i-i+1>\psi R$$

$$<i-i+1>\psi P$$

$$<i-i+1>\psi X$$

$$<1>\psi X$$

$$<1>\psi R$$

$$<1>\psi P$$

$$<1>\psi X$$

$$q\equiv ({}^{-}1△\rho {\xi }_l)[0]\equiv (1△\rho {\xi }_r)[0]$$

$$\rho {\xi }_l{}_{+}{\bullet }_{\times }{\xi }_r\equiv ({}^{-}1\nabla \rho {\xi }_l)++(1\nabla \rho {\xi }_r)$$

$$(<ij>)\psi {\xi }_l({}_{o{p}_0}{\bullet }_{o{p}_1}){\xi }_r\equiv {}_{{}_{o{p}_0}red}(<i>\psi {\xi }_l)o{p}_1((\iota q){}_{++}{\Omega }_{<01>}<j>)\psi {\xi }_r$$

$$q\equiv ({}^{-}1△\rho R)[0]\equiv ({}^{-}1△<2n>)[0]\equiv <n>[0]\equiv n$$

$$(<0>\psi R{)}_{+}{\bullet }_{\times }(<0>\psi R)$$

$$\equiv {}_{{}_{+}red}(\Theta \psi <0>\psi R)\times ((\iota q){}_{++}{\Omega }_{<01>}\Theta )\psi (<0>\psi R)$$

$$\displaystyle\mbox{Here we say, take each scalar, from the left argument $\iota q$},$$

Full text

1. 1 Motivation and Introduction

2. 1.0.1 Intentions of Tutorial

3. 1.0.2 Matrix Multiply and Solvers

4. 2 An Iterative Solver: The Conjugate Gradient (CG)

5. 2.1 Base Case

6. 2.1.1 Pseudocode to MoA: Base Case

7. 2.2 Subsequent Iterations

8. 2.2.1 Pseudocode to MoA: Subsequent Iterations

9. 2.2.2 The Loop Body

10. 2.2.3 Resulting Algorithm

11. 2.2.4 Substituting for $\vec{\alpha}$ and ${\vec{\beta}}$

12. 2.2.5 Transpose

13. 2.2.6 Simplifying

14. 2.3 Psi Reduction

15. 2.3.1 MoA Definition of Inner Product: n-d

16. 2.3.2 Applying the Definition of Inner Product

17. 2.3.3 Another Common Expression

18. 2.3.4 Final Common Section

19. 2.3.5 Putting it All Together

20. 2.4 Complete Iterative Solver, The MoA ONF

21. 2.4.1 Example

22. 3 Translating to Pseudocode

23. 3.1 The CG: ONF to Pseudocode

24. 3.1.1 Implementations from the ONF: Pseudocode

1. 1 Pseudo Code: Base Case
2. 2 MoA: Base Case
3. 3 Pseudo Code: Subsequent Iterations
4. 4 MoA: Subsequent Iterations
5. 5 MoA: Subsequent Iterations - Transpose Removed
6. 6 MoA: Subsequent Iterations - Simplifed Indexing
7. 7 MoA: Subsequent Iterations - Psi Reduced, The ONF
8. 8 Complete Solver: MoA ONF
9. 9 The CG: ONF to Pseudocode

1 Motivation and Introduction

1.0.1 Intentions of Tutorial

It is often difficult to learn new mathematics semantically and syntactically, even when there are similarities in the words and meaning when discussed aloud. The goal of this document is to facilitate learning through explainations and definitions relating our common mathematical knowledge and highlighting what is new. It is meant to be a working document that will evolve based on feedback from target audiences, those mathematically literate in linear and tensor algebra, those that want to learn MoA, Psi Calculus [1], and its uses [2], those that want and need the ability to prove a design, either in hardware or software through the ONF, Operational Normal Form, and those wanting to exploit all resources optimlally, especially when Tensor Algebra, i.e. algorithms foundational to their application, are needed: Knowledge Representation and Modeling, Scientific Computing, Signal Processing, etc.

1.0.2 Matrix Multiply and Solvers

We aim to solve for $\vec{x}$, a vector, in

[TABLE]

where $A$ is symmetric, real, positive-definite, and all variables are conformable for the operaton, i.e. $\vec{x}$, $\vec{b}$, and the number of columns in $A$ have the same number of components.

The matrix-vector product, like other inner products from linear algebra, can be expressed in MoA as ${}_{+}\bullet_{\times}$, denoting point-wise scalar extension of multiplication over the inner dimension, followed by a single application of additive reduction.

We are motivated by a few anamolies in Linear and Tensor Algebra. From a CMU CS professor’s notes on painless-conjugate-gradient [3]: “… A is an n by n matrix, and x and b are vectors, that is, n by 1 matrices. The inner product of two vectors is written $x^{T}y$ and represents the scalar sum $\sum_{i=1}^{n}x_{i}y_{i}$ … In general, expressions that reduce to 1 by 1 matrices are treated as scalars.”

Although we can build symbolic software systems that apply the same equalities as one might do by hand, that implementation is often slow. If we are to have the mathematical prowess of tensor algebra and the computational/communication prowess needed for IoT and the complexity of real time systems with knowledge and reasoning one must realize that it is hard to prove properties of programs if grammars have anamolies. Ideally, if everything was a linear, or multilinear operation, reasoning was possible. We often forget that what can be done in our minds, on paper, etc. mathematically, must now be transformed to an accurate answer done in the most computationally efficient way, given a certain set of resources. What we ignore costs computational resources to analyze and deal with the anamoly. Linear, or multilinear, operations, have the advantage that everything is mathematical without exception. In the above case, at one moment it says a vector is an n by 1, then in the summation it is represented as a vector. They ARE different entities. In a tensor world, a scalar is a 0-dimensional array, a vector is 1-d, and a 1 by 1, is a 2-d array. MoA and Psi Calculus give us this formality while providing the full tensor algebra we are familiar with.

We study the CG because it is the most prominent iterative method for solving sparse systems of linear equations. For dense systems, factor and solve with back substitution is often best. Re complexity, the time spent factoring a dense matrix is roughly equivalent to the time spent solving the system iteratively.

Finallly, we study the CG because, if it is fundamental to other solvers. We aim to paramaterizie the variants, while we assure flexible use, performance, and accuracy.

2 An Iterative Solver: The Conjugate Gradient (CG)

Rewriting the objective with this new notation, we now aim to solve

[TABLE]

where $n$ is a positive integer indicating the size of the problem, and the shapes are denoted as

[TABLE]

and

[TABLE]

2.1 Base Case

The base case populated prior to the first iteration can be expressed by first initializing the guess for $\vec{x}$ as a zero vector of shape $<n>$ (could be set to a better guess if available). Then, we rewrite the pseudocode.

2.1.1 Pseudocode to MoA: Base Case

In the pseudocode we see there are k iterations, implying that a matrix with k rows will hold all recurrences. But, it turns out we only need 2 rows in each matrix used in the recurrence, i.e. the previous time step, and the present time step. We keep the same loop bounds, i.e. $0\leq i<n$, but we change the loop body from $i$ and $i+1$ to $i-i$ and $i+1-i$, i.e. 0 and 1. We then create the following matrices X, P and R, each a 2 by n. Notice how Expression (6) looks now. The Transpose has been eliminated because the transpose of a vector in MoA is the same vector.

2.2 Subsequent Iterations

2.2.1 Pseudocode to MoA: Subsequent Iterations

The iterative solver can be expressed in MoA for a temporal index, $<i>$, where

[TABLE]

effectively replaces the function of $k$ in the pseudocode and allows the procedure to be written using the array equality operator, “$\equiv$”, without requiring the assignment operator, “$:=$”. This can be thought of as unrolling the temporal dimension. This subtlety, is why we are making the change in variable name.

2.2.2 The Loop Body

We rewrite the loop body in MoA syntax by:

Replacing bold vectors with iteration subscripts, with a matrix and an index of that iteration. Vectors that are really vectors are denoted by the same named variable, e.g. x, but as vector symbols: $\mathbf{x}\rightarrow\vec{x}$ 2. 2.

Replacing subscript indexing with the $\psi$ binary function: $\mathbf{r}_{k}\rightarrow<i>\psi{R}$ 3. 3.

Using transpose unary prefix operator $\bigcirc\!\!\!\!\!\backslash\;$: $\mathbf{r}_{k}^{T}\rightarrow\bigcirc\!\!\!\!\!\backslash\;<i>\psi{R}$ 4. 4.

Using inner product binary operator ${}_{+}\bullet_{\times}$, as discussed earlier

2.2.3 Resulting Algorithm

The vectors $\vec{\alpha}$ and $\vec{\beta}$ can be substituted in for each usage. The expanded algorithm is now

2.2.4 Substituting for $\vec{\alpha}$ and ${\vec{\beta}}$

[TABLE]

2.2.5 Transpose

Note that transpose is an unnecessary operation on vectors in MoA (the typical issue of row- and column-vectors encountered in linear algebra is ignored. The reason why, is based on the definition of inner product, which states that,

The shape of the result of the inner product, is the shape of the left argument, exclusive of its last component on the right concatenated to the shape of the right argument, exclusive of its first component. Also, the last component of the left argument, must be the same, as the first component of the right.

So for example, with two vectors, the left argument ${\vec{w}}$, and ${\vec{v}}$, the right argument, both have shape $<5>$, i.e., $\rho\;\vec{w}=<5>$ and $\rho\;\vec{v}=<5>$. Substituting in what was said above:

The last component of the left argument, $(-1)\,\triangle\,(\rho\;\vec{w})=(-1)\,\triangle\,<5>=<5>$, and the first component of the right arguement $1\,\triangle\,(\rho\;\vec{v})=1\,\triangle\,<5>=<5>$, are the same: 2. 2.

Shape of the result: Everything but the last component of the left argument, which is the empty vector, $\Theta$, i. e. $(-1)\,\nabla\,(\rho\;\vec{w})=(-1)\,\nabla\,<5>=\Theta$, concatenated to everything but the first argument of the right argument, also $\Theta$, i. e. $1\,\nabla\,(\rho\;\vec{v})=1\,\nabla\,<5>=\Theta$. Thus, $\Theta+\!\!\!+\Theta$. By definition of concatenation, the concatenation of two empty vectors is $\Theta$, which is the shape of a scalar in MoA, that is, what we know is the result of the inner product of two vectors, or row vectors.

The inner product operator behaves as expected on a pair of vector arguments, so each occurrence of transpose can be eliminated. Also note that, within a single iteration only two index values are used on the LHS of $\psi$: $<i+1>$ and $<i>$. Without loss of generality, we can subtract $i$ from each index to reduce the size of working memory.

2.2.6 Simplifying

2.3 Psi Reduction

2.3.1 MoA Definition of Inner Product: n-d

To paraphrase the definition of inner product, given that $\xi_{l}$ denotes the left argument, and $\xi_{r}$ the right, providing $1\leq\delta\xi_{l}$ and $1\leq\delta\xi_{r}$, i.e. the dimensionality of both arguments must be greater or equal to 1.

[TABLE]

q is a scalar, so we must index the one element vector produced by -1 take of the left shape and 1 take of the right shape. Thus, the shape of the result is defined by

[TABLE]

and for $0\leq i<(\,^{-}1\,\nabla\,\rho\xi_{l})[0]$ and $0\leq j<(\,1\,\nabla\,\rho\xi_{r})[0]$

[TABLE]

2.3.2 Applying the Definition of Inner Product

Before we begin, we observe that many of the pieces of the derivation are the same. Consequently, we pick pieces, reduce to normal form, then put the pieces together to show an ONF for the CG. Applying this definition, we can start to reduce the terms containing relevant inner product applications

[TABLE]

[TABLE]

From the definition of reduction, we know

[TABLE]

We can further simplify expression (33).

[TABLE]

2.3.3 Another Common Expression

The only difference in the following is $<1>\psi R$. Similarly, we can simplify expression  (35).

[TABLE]

2.3.4 Final Common Section

Now, we reduce the remaining inner product term in expressions  (33) and  (34). $\forall\;\;i,j\;\;\mbox{s.t.}\;\;0\leq i<(\rho A)[0]\;\;\mbox{and}\;\;\;0\leq j<(\rho A)[1]$

[TABLE]

2.3.5 Putting it All Together

2.4.1 Example

We consider a linear system $A_{+}\bullet_{\times}\vec{x}\equiv\vec{b}$ given by

[TABLE]

This system differs from the linear algebra formulation by removing details about row and column vectors that are unnecessary for finding a solution. While we could start with any guess without loss of generality, we will consider

[TABLE]

The remaining initial values are then

[TABLE]

Since $0\leq k<n$, we can expand easily for all values $k$

[TABLE]

So our initial condition is

[TABLE]

We now evaluate the next approximate solution by solving for $<1>\psi X$ for all valid $k$ and find

[TABLE]

Substituting in the calculated initial condition, we have

[TABLE]

which evaluates to

[TABLE]

We then compute the next residual vector $<1>\psi R$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The new residual vector is now $<1>\psi R\equiv<\,^{-}0.2810\;0.7492>$. Assuming this does not meet our convergence criterion, we will continue to calculate $<1>\psi P$.

[TABLE]

So, $<1>\psi P\equiv<\,^{-}0.3512\;0.7229>$. Now we can iterate with the new values of $X$, $R$, and $P$ until a satisfactory solution is obtained.

3 Translating to Pseudocode

Note that we have switched back to using an assignment operator “:=” to mimic execution behavior. Recall, that this definition can be readily translated to a C-like language by replacing each raveled array with a pointer to an address in memory, denoted by X, R. and P, and using the square brackets as the C array dereference operator.

3.1 The CG: ONF to Pseudocode

3.1.1 Implementations from the ONF: Pseudocode

We now have a design for hardware or software.