Clojure Linear Algebra Refresher (1) - Vector Spaces

To work with vectors in a vector space \(R^n\), I use Neanderthal vectors, and specify the dimension \(n\) as a part of the construction of the vector. Sometimes the general vctr function is appropriate, and sometimes it is easier to work with specialized functions such as dv and sv (for native vectors on the CPU), clv (OpenCL vectors on GPU or CPU), or cuv (CUDA GPU vectors). This tutorial will use native, double floating-point precision vectors on the CPU (dv).

To get access to these constructors, I'll require the namespace that defines them. I'll also require the core namespace, where most functions that operate on vectors and matrices are located.

(require '[uncomplicate.neanderthal
           [native :refer :all]
           [core :refer :all]])

For example, vector v1 is a vector in \(R^4\), while v2 is a vector in \(R^{22}\).

(def v1 (dv -1 2 5.2 0))
(def v2 (dv (range 22)))

If i print these vectors out in the REPL, I can find more details about them, including their components (elements, entries):

v1

#RealBlockVector[double, n:4, offset: 0, stride:1]
[  -1.00    2.00    5.20    0.00 ]

v2

#RealBlockVector[double, n:22, offset: 0, stride:1]
[   0.00    1.00    2.00    ⋯     20.00   21.00 ]

Dot product of two compatible vectors

\(\mathbf{u}\cdot\mathbf{v} = u_1 v_1 + u_2 v_2 + \cdots + u_n v_n\)

can be computed using the function with a predictably boring name dot:

(let [u (dv 1 -2 4)
      v (dv 3 0 2)]
  (dot u v))

11.0

Norm (aka length, or magnitude, or L2 norm) most often refers to Euclidean distance between two vectors (although other norms can be used):

\(\lVert \mathbf{u} \lVert = \sqrt{u_1^2 + u_2^2 + \cdots + u_n^2}\)

In Clojure, we use the nrm2 function to compute this norm:

(nrm2 (dv 1 3 5))

5.916079783099616

A unit vector is a vector whose norm is 1. Given a vector v, we can construct a unit vector u in the same direction:

\(\mathbf{u} = \frac{1}{\lVert \mathbf{v} \lVert} \mathbf{v}\)

or, in Clojure code:

(let [v (dv 1 3 5)]
  (scal (/ (nrm2 v)) v))

#RealBlockVector[double, n:3, offset: 0, stride:1]
[   0.17    0.51    0.85 ]

This process is called normalizing the vector.

It is important to note that real vectors are not the only "kind" of vector spaces there is. Anything that has operations of addition and scalar multiplication that have certain properties (see the textbooks for details) is a vector space. Some well known vector spaces are real vectors (\(R^n\)), matrices (\(M^{mn}\)), complex vectors (\(C^n\)), and, of course - functions. However, when we do numerical computing, we usually work with real or complex vectors and matrices, so in the following discussion I'll refer only to the parts of the textbook that deal with those. For the abstract theoretical parts, you won't use Clojure anyway, but your trusty pen and paper.

Neanderthal currently supports only real vectors and matrices; complex vectors and matrices will be supported in some of the upcoming versions.

I have already shown you how to use vectors. With matrices, it's similar, just with more options. The first difference is that we can use different optimized implementations of matrices, to exploit knowledge of their structure. Therefore, there are dense matrices (GE), dense triangular matrices (TR), dense symmetrical matrices (SY), various banded storages for symmetric and triangular matrices, sparse matrices (Sp), etc. Some of those are already implemented and supported in Neanderthal at the time of this writing, some are waiting their turn to be implemented (and may be at the time you're reading this - check Neanderthal API docs out).

Anyway, let's construct a matrix and compute its norm:

(def m (dge 2 3 (range 6)))
m

#'user2/m#RealGEMatrix[double, mxn:2x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓
   →       0.00    2.00    4.00
   →       1.00    3.00    5.00
   ┗                               ┛

(nrm2 m)

7.416198487095662

Euclidean norm applied to matrices is known as Frobenius norm (FYI).

Most of the functions that can be applied to vectors, can also be applied to matrices. Furthermore, there are many more functions in linear algebra that work exclusively with matrices. Neanderthal API docs list what is available in Clojure at the moment when you read this.

The following two sections, Subspaces and Linear Combinations are theoretical. We skip them, since I can not show you anything in Clojure code related to this.

Lots of theory here, but when it comes to computation, I can note that two vectors are linearly dependent if they are on the same line, implying that the angle between them is 0, meaning the cosine is 1. If cosine is not 1, they are independent. There is a catch, though: floating-point computations are not absolutely precise. For very large vectors, sometime there will be rounding error, so you might get something like 0.99999999 or 1.000001. Equality comparison of floating-point numbers is a tricky business. Keep that in mind, and know that even there Neanderthal can help you, since it offers functions that take some imprecision margin into account. See the Neanderthal API docs and look up for functions such as f= in the math namespace.

So, in Clojure:

(let [x (dv 2 -1 3)
      y (dv 4 -2 6)]
  (/ (dot x y) (nrm2 x) (nrm2 y)))

1.0

These two vectors are linearly dependent indeed.

But, what if we have more than two vectors, and are tasked with ckecking whether their set is dependent or not? See example 1 from the page 170 of the book I recommended. There is a set of vectors \({x = (1,2,0), y = (0,1,-1), z = (1,1,2)}\), and we need to find whether there are scalars \(c_1, c_2, \text{and } c_3\) such that \(c_1\mathbf{x} + c_2\mathbf{y} + c_3\mathbf{z} = 0\). This leads to a system of linear equations that, when solved, could have an unique solution, in which case the vectors are linearly independent, or to more than one solution, in which case the vectors are linearly dependent.

Luckily for us, Neanderhal offers functions for solving systems of linear equations. Require the linalg namespace to reach for them:

(require '[uncomplicate.neanderthal.linalg :refer :all])

Now, we can set up and solve the system that will find the answer for us:

(let [cs (dge 3 1)]
  (sv! (dge 3 3 [1 2 0
                 0 1 -1
                 1 1 2])
       cs)
  cs)

#RealGEMatrix[double, mxn:3x1, layout:column, offset:0]
   ▥       ↓       ┓
   →       0.00
   →      -0.00
   →      -0.00
   ┗               ┛

Just like in the textbook, the solution is: \(c_1=0, c_2=0, c_3=0\). If we tried with obviously dependent vectors, we would have got an exception complaining that the solution could not be computed. Not bad…

Some more theory here that does not relate directly to Clojure code.

Note that the dimension of the space of a vector could be inquired with the (dim x) function. The dimension of the space of the matrix is \(m \times n\), so (* (mrows a) (ncols a)). All the spaces we deal with here have finite dimension, while we leave infinite dimension vectors to pen and paper adventures.

Rank of a matrix is useful since it relates matrices to vectors, and can tell us some useful properties of a matrix at hand. See the textbook (or wikipedia) for the definition and use. What is the question right now is how do we compute it? It is not easy at all, but, fortunately, I can use Neanderthal. When we inquire about the rank, we expect a natural number as a result.

We can get rank by doing singular value decomposition by calling function svd, and inspecting one of its results, namely the sorted vector of singular values s. The number of singular values in :sigma that are greater than 0 is the rank of the matrix. As before, take into account that floating point numbers should be compared with equality carefully.

Here's example 3 (page 185 from the textbook) done in Clojure:

(let [a (dge 4 4 [1 0 0 0
                  2 0 0 0
                  0 1 0 0
                  0 0 1 0])]
  (svd a))

#uncomplicate.neanderthal.internal.common.SVDecomposition{:sigma #RealDiagonalMatrix[double, type:gd mxn:4x4, offset:0]
   ▧                                       ┓
   ↘       2.24    1.00    1.00    0.00
   ┗                                       ┛
, :u nil, :vt nil, :master true}

After inspecting vector :sigma for non-zero values (visually, or by writing a small function that checks its elements) I am a little smarter than before since I know now that \(rank(A) = 3\).

Seems so cumbersome. Why isn't there a function rank in Neanderthal that does this for me? Well, the problem is that I rarely need only rank. Often I also need those other results from svd! or similar functions. If rank was so easy to access, it would have promoted wasteful computations of other results that go with it. If you really need rank it is easy to write that function as part of your personal toolbox. Maybe I'll include such functions in Neanderthal later. Who knows. What is important is that the cake is available right now, so you can easily add your own cherries on top of it.

This is yet another theoretical section. From the computational point of view, what could be interesting, is computing projection of a vector onto another vector or a subspace.

\(proj_u v = \frac{v \cdot u}{u \cdot u} u\)

Basically, it boils down to be able to compute dot product and then scale vector u.

(let [v (dv 6 7)
      u (dv 1 4)]
  (scal (/ (dot u v) (dot u u)) u))

#RealBlockVector[double, n:2, offset: 0, stride:1]
[   2.00    8.00 ]

It's similar for subspaces, so that might be the obligatory trivial exercise that is being left to the reader.

That's all folks… For this post, at least. I've covered most (or perhaps all) numerical stuff that is typically covered in a Vector Spaces chapter of a typical textbook. I am not a mathematician, so if you are, please forgive me for not being precise enough, and if I've made errors (I'm sure I have), please send me corrections.

I hope that this post helped you started with doing linear algebra computations in Clojure. You might also want to continue reading the next article: Eigenvalues and Eigenvectors.