Linear Algebra

This tutorial covers vector and matrix representations and concepts, including vector slices, matrix blocks, and least-squares solutions.

Amino performs linear algebra using the standard, high-performance-computing BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra PACKage) libraries. Multiple, optimized implementations of BLAS/LAPACK are available such as ATLAS, OpenBLAS, and the Intel Math Kernel Library (MKL).

Vectors

Representation

We represent mathematical vectors using a length count, a data pointer, and an "increment" between successive elements. For example, consider the following vector:

\[ \mathbf{x} = \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{len:} & 6 \\ \texttt{data:} & \begin{array}{|c|c|} \hline 1 & 2 & 3 & 4 & 5 & 6\\ \hline \end{array} \\ \texttt{inc:} & 1 \\ \end{cases} \]

We can "slice" every other element of \(\mathbf{x}\) by using the same data pointer and an increment of 2. The length is now 3.

\[ \mathbf{x}_{0:6:2} = \begin{bmatrix} 1 \\ 3 \\ 5 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{len:} & 3 \\ \texttt{data:} & \begin{array}{|c|c|} \hline 1 & 2 & 3 & 4 & 5 & 6\\ \hline \end{array} \\ \texttt{inc:} & 2 \\ \end{cases} \]

The combination of length, data pointer, and increment value items enables flexible access to parts of other vectors and matrices, such as "slices" of vectors and rows, columns, or diagonals of matrices. The following C struct defines the vector for doubles:

struct aa_dvec {
    size_t len;
    double *data;
    size_t inc;
}

Example Code

The following code illustrates the vector representations and slices.

#!/usr/bin/env python3
 
from amino import DVec
 
# Create a Vector
x = DVec([1, 2, 3, 4])
print("x:     %s" % x)
print("x.len: %d" % len(x))
print("x.inc: %d" % x.inc)
 
# Common Arithmetic Operators
print("-x:   %s" % (-x))
print("2*x:  %s" % (2 * x))
print("x+x:  %s" % (x + x))
print("x+1:  %s" % (x + 1))
print("x/2:  %s" % (x / 2))
 
# Slice the vector
s = x[0:4:2]
print("s: %s" % s)
print("s.len: %d" % len(s))
print("s.inc: %d" % s.inc)
 
# Increment the Slice
s += 1
print("s: %s" % s)
print("x: %s" % x)

Matrices

Representation

The memory layout for matrices stores elements column-by-column, i.e., column-major format. For example, a 3x3 matrix would be stored in memory as follows:

\[ \mathbf{A} = \begin{bmatrix} x_0 & y_0 & z_0 \\ x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{rows:} & 3 \\ \texttt{cols:} & 3 \\ \texttt{data:} & \begin{array}{|c|c|} \hline x_0 & x_1 & x_2 & y_0 & y_1 & y_2 & z_0 & z_1 & z_2 \\ \hline \end{array} \end{cases} \]

Thus, to represent the entire matrix, we need the row count, the column count, and a pointer to the data. To represent a block or submatrix, we need one additional element: the "leading dimension," i.e., the number of rows in the containing matrix. Consider the block consisting of the first two rows and columns of \(\mathbf{A}\). We use the same data pointer, but the distance between successive columns of this block is NOT the row count of the block, but the row count of overall matrix \(\mathbf{A}\):

\[ \mathbf{A}_{0:2,0:2} = \begin{bmatrix} x_0 & y_0 \\ x_1 & y_1 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{rows:} & 2 \\ \texttt{cols:} & 2 \\ \texttt{data:} & \underbrace{ \begin{array}{|c|c|} \hline x_0 & \strut x_1 & x_2 \\ \hline \end{array} }_{\texttt{ld=3}} \!\begin{array}{|c|c|} \hline y_0 & \strut y_1 & y_2 \\ \hline \end{array} \!\!\!\!\ \begin{array}{|c|c|} \hline z_0 & \strut z_1 & z_2 \\ \hline \end{array} \\ \texttt{ld:} & 3 \end{cases} \]

Thus, our matrix representation consists of the row count, column count, data pointer, and leading dimension. Note that for an entire matrix, the leading dimension will equal the row count.

struct aa_dmat {
    size_t rows;
    size_t cols;
    double *data;
    size_t ld;
}

We can use the previously introduced vector format to view rows, columns, and diagonals of our matrix.

A column vector has length equal to the number of rows, a data pointer to the first element of the column, and a increment of 1:

\[ \texttt{col}(\mathbf{A},0) = \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{len:} & 3 \\ \texttt{data:} & \begin{array}{|c|c|} \hline x_0 & x_1 & x_2 \\ \hline \end{array} \\ \texttt{inc:} & 1 \end{cases} \]

\[ \texttt{col}(\mathbf{A},1) = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{len:} & 3 \\ \texttt{data:} & \begin{array}{|c|c|} \hline y_0 & y_1 & y_2 \\ \hline \end{array} \\ \texttt{inc:} & 1 \end{cases} \]

A row vector has length equal to the number of columns, a data pointer to the first element in the row, and an increment equal to the matrix leading dimension.

\[ \texttt{row}(\mathbf{A},0) = \begin{bmatrix} x_0 & y_0 & z_0 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{len:} & 3 \\ \texttt{data:} & \begin{array}{|c|c|} \hline x_0 & x_1 & x_2 & y_0 & y_1 & y_2 & z_0 \\ \hline \end{array} \\ \texttt{inc:} & 3 \end{cases} \]

\[ \texttt{row}(\mathbf{A},1) = \begin{bmatrix} x_1 & y_1 & z_1 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{len:} & 3 \\ \texttt{data:} & \begin{array}{|c|c|} \hline x_1 & x_2 & y_0 & y_1 & y_2 & z_0 & z_1 \\ \hline \end{array} \\ \texttt{inc:} & 3 \end{cases} \]

The matrix diagonal vector has length equal to the minimum of rows and columns, data pointer to the matrix data, and increment equal to one plus the leading dimension:

\[ \texttt{diag}(\mathbf{A}) = \begin{bmatrix} x_0 \\ y_1 \\ z_2 \end{bmatrix} \qquad \leadsto \qquad \begin{cases} \texttt{len:} & 3 \\ \texttt{data:} & \begin{array}{|c|c|} \hline x_0 & x_1 & x_2 & y_0 & y_1 & y_2 & z_0 & z_1 & z_2 \\ \hline \end{array} \\ \texttt{inc:} & 4 \end{cases} \]

Example Code

#!/usr/bin/env python3
 
from amino import DMat
 
 
# Helpers to print(headings)
def h1(name):
    """Print a level 1 heading"""
    print("")
    print("{:^16}".format(name))
    print("{:=^16}".format(''))
 
 
def h2(name):
    """Print a level 2 heading"""
    print("")
    print("{:^16}".format(name))
    print("{:-^16}".format(''))
 
 
# Create a 2x3 matrix
h1("A")
A = DMat.row_matrix([[1, 2, 3], [4, 5, 6]])
print(A)
print("A.rows: %d" % A.rows)
print("A.cols: %d" % A.cols)
print("A.ld:   %d" % A.ld)
for k in range(6):
    print("A._data[%d] = %f" % (k, A._data[k]))
 
# Create a 3x2 matrix
h1("B")
B = DMat.col_matrix([[1, 2, 3], [4, 5, 6]])
print(B)
print("B.rows: %d" % B.rows)
print("B.cols: %d" % B.cols)
print("B.ld:   %d" % B.ld)
for k in range(6):
    print("B._data[%d] = %f" % (k, A._data[k]))
 
# Example matrix arithmetic
h1("Arithmetic")
 
h2("A*2")
print(A * 2)
 
h2("A/2")
print(A / 2)
 
h2("A+A")
print(A + A)
 
h2("A+1")
print(A + 1)
 
h2("A-1")
print(A - 1)
 
h2("A^T")
print(A.transpose())
 
h2("A*[1,2,3]^T")
print(A * [1, 2, 3])
 
h2("A*B")
print(A * B)
 
# Vector views of matrix
h1("Matrix Vector Views")
 
h2("A")
print(A)
 
h2("A.row_vec(0)")
print(A.row_vec(0))
 
h2("A.col_vec(1)")
print(A.col_vec(1))
 
h2("A.diag_vec()")
d = A.diag_vec()
print(d)
 
h2("A.diag_vec() += 10")
d += 10
print(d)
 
h2("A")
print(A)
 
# Block views of matrix
h1("Matrix Block Views")
 
h2("A")
print(A)
 
h2("As = A[0:2,0:2]")
As = A[0:2, 0:2]
print(As)
print("As.rows: %d" % As.rows)
print("As.cols: %d" % As.cols)
print("As.ld:   %d" % As.ld)
 
h2("As *= 2")
As *= 2
print(As)
 
h2("A")
print(A)
 
h2("B")
print(B)
 
h2("Bs = B[0:2,0:2]")
Bs = B[0:2, 0:2]
print(Bs)
print("Bs.rows: %d" % Bs.rows)
print("Bs.cols: %d" % Bs.cols)
print("Bs.ld:   %d" % Bs.ld)
 
h2("Bs *= 2")
Bs *= 2
print(Bs)
 
h2("B")
print(B)

Linear Least Squares

This sections discusses how to solve systems of linear equations using the matrix inverse and pseudoinverse.

(An informed reader may note that explicitly computing the (pseudo)inverse is not the most efficient method for solving linear equations. However, obtaining an explicit pseudoinverse provides useful capabilities when we are solving for a robot's motion, so we will present the explicit computation here.)

Systems of Linear Equations

We can represent and efficiently solve systems of linear equations in terms of matrices. Consider the example below with two linear equations and two unknowns. We represent the linear coefficients of each unknown as a matrix.

\[ \begin{array} 1 x_0 & + & 2 x_1 & = 5 \\ 3 x_0 & + & 4 x_1 & = 6 \end{array} \quad\quad \leadsto \quad\quad \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \end{bmatrix} = \begin{bmatrix} 5 \\ 6 \end{bmatrix} \]

Then, we can solve the system of equations by multiplying each side by the inverse of the coefficient matrix.

\[ \begin{bmatrix} x_0 \\ x_1 \end{bmatrix} = \left( \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}^{-1} \right) \begin{bmatrix} 5 \\ 6 \end{bmatrix} = \begin{bmatrix} -2 & 1 \\ 1.5 & -.5 \end{bmatrix} \begin{bmatrix} 5 \\ 6 \end{bmatrix} = \begin{bmatrix} 4 \\ 4.5 \end{bmatrix} \]

Over and under-determined Systems

If we have more unknowns than equations, there may be infinitely many solutions to the system. Similarly, more equations than unknowns may result in a system with no solutions. We can still construct a coefficient matrix, but it will now be nonsquare:

In the underdetermined case, we have a "fat" matrix (more rows than columns):

\[ \begin{array} 1 x_0 & + & 2 x_1 & 3 x_2 & = & 7 \\ 4 x_0 & + & 5 x_1 & 6 x_2 & = & 8 \end{array} \quad\quad \leadsto \quad\quad \begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 7 \\ 8 \end{bmatrix} \]

In the overdetermined case, we have a "thin" matrix (more columns than rows:

\[ \begin{array} 1 x_0 & + & 2 x_1 \\ 3 x_0 & + & 4 x_1 \\ 5 x_0 & + & 6 x_1 \end{array} \quad\quad \leadsto \quad\quad \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix} \begin{bmatrix} y_0 \\ y_1 \end{bmatrix} = \begin{bmatrix} 7 \\ 8 \\ 9 \end{bmatrix} \]

Since we have a nonsquare matrix, we cannot compute it's inverse. However, we can use the pseudoinverse.

\[ \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \left( \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix}^{+} \right) \begin{bmatrix} 7 \\ 8 \end{bmatrix} = \begin{bmatrix} -9.44 & 0.444 \\ -0.111 & 0.111 \\ 0.722 & -0.222 \end{bmatrix} \begin{bmatrix} 7 \\ 8 \end{bmatrix} = \begin{bmatrix} -3.06 \\ 0.111 \\ 3.28 \end{bmatrix} \]

\[ \begin{bmatrix} y_0 \\ y_1 \end{bmatrix} \approx \left( \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix}^{+} \right) \begin{bmatrix} 7 \\ 8 \\ 9 \end{bmatrix} = \begin{bmatrix} -1.33 & -0.333 & 0.667 \\ 1.08 & 0.333 & -0.417 \end{bmatrix} \begin{bmatrix} 7 \\ 8 \\ 9 \end{bmatrix} = \begin{bmatrix} -6 \\ -6.5 \end{bmatrix} \]

The pseudoinverse finds the least squares solution.

In the underdetermined case, we are minimizing the sum of squares (i.e., the Euclidean norm) of the vector \(\mathbf{x}\):

• Given: Underdetermined system \(\mathbf{A}\mathbf{x} = \mathbf{b}\), i.e., \(\mathbf{A}\) has more columns (unknowns) than rows (equations).
• Find: Vector \(\mathbf{x}\), such that:
  • minimize: \(\Vert\mathbf{x}\Vert\)
  • subject to: \(\mathbf{A} \mathbf{x} = \mathbf{b}\)
• Solution: \( \mathbf{x} = \mathbf{A}^+ \mathbf{b} \)

In the overdetermined case, we are minimizing the sum of squares of the error:

• Given: Overdetermined system \(\mathbf{A}\mathbf{x} = \mathbf{b}\), , i.e., \(\mathbf{A}\) has more rows (equations) than columns (unknowns).
• Find: Vector \(\mathbf{z}\), such that:
  • minimize: \(\Vert\mathbf{A} \mathbf{z} - \mathbf{b}\Vert\)
• Solution: \( \mathbf{x} = \mathbf{A}^+ \mathbf{b} \)

Example Code

#!/usr/bin/env python3
 
from amino import DMat
 
 
def h1(name):
    """Print a level 1 heading"""
    print("")
    print("{:^16}".format(name))
    print("{:=^16}".format(''))
 
 
def h2(name):
    """Print a level 2 heading"""
    print("")
    print("{:^16}".format(name))
    print("{:-^16}".format(''))
 
 
def result(A, x, b):
    h2("A")
    print(A)
    h2("b")
    print(b)
    h2("x")
    print(x)
 
 
h1("Fully-Determined: A*x = b")
A = DMat.row_matrix([[1, 2], [3, 4]])
b = [5, 6]
x = A.inv() * b
result(A, x, b)
 
h1("Over-Determined: A*x = b")
A = DMat.row_matrix([[1, 2], [3, 4], [5, 6]])
b = [7, 8, 9]
x = A.pinv() * b
result(A, x, b)
 
h1("Under-Determined: A*x = b")
A = DMat.row_matrix([[1, 2, 3], [4, 5, 6]])
b = [7, 8]
x = A.pinv() * b
result(A, x, b)

See Also

Python

• amino.mat.DVec
• amino.mat.DMat

C

• mat.h