The goal of this post is to be able to solve equations like this one:
\[\left(\begin{array}{cccc}
0.18& 0.60& 0.57& 0.96
0.41& 0.24& 0.99& 0.58 \
0.14& 0.30& 0.97& 0.66 \
0.51& 0.13& 0.19& 0.85 \end{array} \right )
\left ( \begin{array}{c} x_0 \ x_1 \ x_2 \ x_3 \end{array} \right )
= \left ( \begin{array}{c} 1.0 \ 2.0 \ 3.0 \ 4.0 \end{array} \right ) \]
The answer is given as
\[x = \left ( \begin{array}{c} -4.05205 \ -12.6056 \ 1.66091 \ 8.69377 \end{array} \right ) \]
The syntax we want to use is shown below, and we want it to return a vector containing the solution:
(let ((A [[0.18 0.60 0.57 0.96]
[0.41 0.24 0.99 0.58]
[0.14 0.30 0.97 0.66]
[0.51 0.13 0.19 0.85]])
(b [1.0 2.0 3.0 4.0]))
(gsl-linalg-LU-solve A b))Rather than put all the code in here like I have for the past several posts, I started a git repo at https://github.com/jkitchin/emacs-modules that contains this code.
The module for this post can be found here: https://github.com/jkitchin/emacs-modules/blob/master/gsl-linalg.c. There are a few notable features in it. First, I started writing/collecting some helper functions to make these modules simpler to write. For example, look how nice this looks to declare the functions and provide the feature.
int emacs_module_init(struct emacs_runtime *ert)
{
emacs_env *env = ert->get_environment(ert);
DEFUN("gsl-linalg-LU-solve", Fgsl_linalg_LU_solve, 2, 2,
"(gsl-linalg-LU-solve A b).\n" \
"Solve A x = b for x.\n" \
"Returns a vector containing the solution x.",
NULL);
provide(env, "gsl-linalg");
return 0;
}The DEFUN and provide function are defined in https://github.com/jkitchin/emacs-modules/blob/master/emacs-module-helpers.c.
Within the module itself, we have to loop over the inputs to create the arrays that GSL wants to solve this problem. Second, after the solution is obtained, we have to build up a vector to return. The solution is in a gsl_vector, and we need to create an array of emacs_value elements containing each element of the gsl_vector as a float, and then create a vector to return to emacs. I use vectors here because it was easy to get the size of the b vector, which is also related to the size of the A matrix.
The repo has a Makefile in it, so we can build this module with:
make gsl-linalg.soOnce it is compiled, we load it like this. In this post we are in the emacs-modules directory where the gsl-linalg.so library is, and it is not on my load-path, so I add it here.
(add-to-list 'load-path (expand-file-name "."))
(require 'gsl-linalg)Here is one function in the module:
(describe-function 'gsl-linalg-LU-solve)Now, we can solve linear equations like this:
(gsl-linalg-LU-solve
[[0.18 0.60 0.57 0.96]
[0.41 0.24 0.99 0.58]
[0.14 0.30 0.97 0.66]
[0.51 0.13 0.19 0.85]]
[1.0 2.0 3.0 4.0])We have a limited ability to confirm this answer. I have written a function that uses blas for multiplication of 2d vectors. You can see from this:
(gsl-blas-dgemm [[0.18 0.60 0.57 0.96]
[0.41 0.24 0.99 0.58]
[0.14 0.30 0.97 0.66]
[0.51 0.13 0.19 0.85]]
[[-4.052050229573973]
[-12.605611395906903]
[1.6609116267088417]
[8.693766928795227]])That within float that indeed
The main limitation of this module at the moment is that you have to use vectors; you cannot put in a list of numbers. It is possible to make it take lists and vectors, but for now I am leaving it at vectors. Also, it only produces solutions of float numbers (not integers).
The module does not handle 1d vectors well,, e.g. in gsl-linalg-LU-solve example, the right hand side is implied to be a column vector, and we don’t have the array broadcasting features of Python yet. Those are doable things for some future day perhaps. For now I am happy to have figured out how to handle arrays!
In a previous post I implemented a Newton solver in elisp to solve some problems numerically. Today, we continue the dynamic module studies and implement a bracketed root solver from the GNU Scientific Library (https://www.gnu.org/software/gsl/doc/html/roots.html#examples). I will implement the Brent solver here, which is a bracketed root finder: You specify a bracket that contains the root, and the software automatically finds it. It isn’t my favorite solver, I prefer a single initial guess, but this one was easier to implement for now, and the root polishing algorithms in GSL seem to all require the function and derivatives, which I did not want to get into now.
The elisp signature I want to solve an equation
(gsl-root-fsolver-brent f xlo xhi &optional params epsabs epsrel)So, here it is in action.
(add-to-list 'load-path (expand-file-name "."))
(require 'gsl-roots)Here is a simple equation
(gsl-root-fsolver-brent (lambda (x params) (- (* x x) 5)) 0.0 5.0)For comparison:
(sqrt 5)These differ in about the 5th decimal place. If we lower the relative error (the default is only 1e-3), we get quantitative agreement with the analytical solution.
(gsl-root-fsolver-brent (lambda (x params) (- (* x x) 5)) 0.0 5.0 nil nil 1e-6)This module
#include "emacs-module.h"
#include <stdlib.h>
/* Declare mandatory GPL symbol. */
int plugin_is_GPL_compatible;
/* Bind NAME to FUN. */
static void bind_function (emacs_env *env, const char *name, emacs_value Sfun)
{
/* Set the function cell of the symbol named NAME to SFUN using
the 'fset' function. */
/* Convert the strings to symbols by interning them */
emacs_value Qfset = env->intern (env, "fset");
emacs_value Qsym = env->intern (env, name);
/* Prepare the arguments array */
emacs_value args[] = { Qsym, Sfun };
/* Make the call (2 == nb of arguments) */
env->funcall (env, Qfset, 2, args);
}
/* Provide FEATURE to Emacs. */
static void
provide (emacs_env *env, const char *feature)
{
/* call 'provide' with FEATURE converted to a symbol */
emacs_value Qfeat = env->intern (env, feature);
emacs_value Qprovide = env->intern (env, "provide");
emacs_value args[] = { Qfeat };
env->funcall (env, Qprovide, 1, args);
}
/* */
// This just returns the argument, works for a list.
static emacs_value Ff1 (emacs_env *env, int nargs, emacs_value args[], void *data)
{
return args[0];
}
// get first element of a vector
static emacs_value Ff2 (emacs_env *env, int nargs, emacs_value args[], void *data)
{
return env->vec_get(env, args[0], 0);
}
// This just returns a vector of integers!!!
static emacs_value Ff3 (emacs_env *env, int nargs, emacs_value args[], void *data)
{
int len = 2;
emacs_value *array = malloc(sizeof(emacs_value) * len);
array[0] = env->make_integer(env, 2);
array[1] = env->make_integer(env, 4);
emacs_value Fvector = env->intern(env, "vector");
emacs_value vec = env->funcall(env, Fvector, len, array);
free(array);
return vec;
}
// return vector * n
static emacs_value Ff4 (emacs_env *env, int nargs, emacs_value args[], void *data)
{
emacs_value input = args[0];
double N = env->extract_float(env, args[1]);
int len = env->vec_size (env, input);
emacs_value *array = malloc(sizeof(emacs_value) * len);
// multiply each value by N
for (ptrdiff_t i = 0; i < len; i++)
{
array[i] = env->make_float(env,
N * env->extract_float(env,
env->vec_get (env, input, i)));
}
// If you change this to list, you get a list instead!
emacs_value Fvector = env->intern(env, "vector");
emacs_value vec = env->funcall(env, Fvector, len, array);
free(array);
return vec;
}
// return 2nd element of vector
static emacs_value Ff5 (emacs_env *env, int nargs, emacs_value args[], void *data)
{
emacs_value vec = args[0];
return env->vec_get (env, vec, 1);
}
// get second value of second vector
static emacs_value Ff6 (emacs_env *env, int nargs, emacs_value args[], void *data)
{
emacs_value vec = args[0];
emacs_value v2 = env->vec_get (env, vec, 1);
return env->vec_get (env, v2, 1);
}
// index a list
static emacs_value Ff7 (emacs_env *env, int nargs, emacs_value args[], void *data)
{
emacs_value nth = env->intern(env, "nth");
return env->funcall (env, nth, 2, args);
}
int emacs_module_init (struct emacs_runtime *ert)
{
emacs_env *env = ert->get_environment (ert);
#define DEFUN(lsym, csym, amin, amax, doc, data) \
bind_function (env, lsym, \
env->make_function (env, amin, amax, csym, doc, data))
DEFUN("f1", Ff1, 1, 1, NULL, NULL);
DEFUN("f2", Ff2, 1, 1, NULL, NULL);
DEFUN("f3", Ff3, 0, 0, NULL, NULL);
DEFUN("f4", Ff4, 2, 2, NULL, NULL);
DEFUN("f5", Ff5, 1, 1, NULL, NULL);
DEFUN("f6", Ff6, 1, 1, NULL, NULL);
DEFUN("f7", Ff7, 2, 2, NULL, NULL);
provide (env, "mod-vector");
/* loaded successfully */
return 0;
}rm -f mod-vector.so mod-vector.o
gcc -Wall -I/usr/local/include -fPIC -c mod-vector.c
gcc -shared -L/usr/local/include -lgsl -o mod-vector.so mod-vector.o(org-babel-tangle)In a previous post I integrated some linear algebra into Emacs using the GNU Scientific library and a dynamic module. In this post, I use a similar approach that uses the Intel MKL library in conjunction with some helper elisp functions to mimic the array broadcasting features in Numpy. I thought this might be easier and lead to at least a complementary set of functionalities.
Note: I had to follow the directions here to disable some security feature on my Mac so that it would use the MKL libraries. Thanks Apple.
It is convenient to use vectors for the representation of arrays in Emacs because there are nice functions in the emacs-module.h for accessing vector properties. Also vectors sound closer to an array than a list. So what about array broadcasting, e.g. the way numpy lets you multiply a 2d array with a 1d array? If you multiply two arrays with size (m1, n1) * (m2, n2), it is required that the number of columns in the first array (n1) be equal to the number of rows in the second one (m2), and the resulting size of the array product will be (m1, n2). What should happen though when we have 1d array? This is neither a row or column vector itself, but we can treat as either one if we choose too. For example the vector [1 2 3] can be thought of as an array with the shape (1, 3), e.g. a single row with three columns, or (3, 1), i.e. three rows in a single column. We will build this capability into the module for convenience.
I still find it moderately tedious to write c functions that take emacs arguments, transform them to c arguments, do some c computations, and convert the results back to emacs values. So, we only implement one c function for this that multiplies two 2d arrays together using the cblas_dgemm routine in the MKL library. Then, we will create a complementary elisp library that will provide some additional functionality to get the shapes of vector arrays, dimensions, and allow us to multiply 1d and 2d vectors together the same way Numpy does array broadcasting.
The dynamic module code is listed in The module code. The elisp code is listed in Elisp helper functions. In the following sections we just demonstrate how to use the results.
I found it convenient to do array shape and dimension analysis in elisp before sending arrays to the dynamic module. The shape of an array is just the number of elements in each dimension. Here we look at a 2× 3 array.
(vector-shape [[1 2 3]
[3 4 5]])You see it returns a vector showing two rows and three columns. There are two convenience commands to get the number of rows (vector-nrows) and columns (vector-ncols). Here is one of them.
(vector-ncols [[1 2 3]
[3 4 5]])The main problem we want to calculate is the matrix multiplication
This is a simple dot-product that is actually calculated in elisp.
(matrix-multiply [1 1 1] [2 2 2])✓
Note we get a float. That is because we initialize the sum with 0.0 to be consistent with all the other cases which are done with floats. dgemm is a double routine in MKL, which means it should return floats. Internally in the module, we cast all numbers as doubles for the multiplication.
This is a matrix multiplication that is typically like
\[ \left[\begin{array}{cc}
1 & 2
3 & 4 \end{array}\right]
\left [ \begin{array}{c}
1 \ 1 \end{array}\right] = \left[\begin{array}{c}3\7\end{array}\right]\]
(let ((A [[1 2]
[3 4]])
(b [1 1]))
(matrix-multiply A b))✓
This case is
\[\left[\begin{array}{cc}1 & 1\end{array}\right] \left[\begin{array}{cc} 1 & 2\ 3 & 4\end{array}\right] = \left[\begin{array}{cc} 4 & 6 \end{array}\right ]\]
(matrix-multiply [1 1]
[[1 2]
[3 4]])✓
As with the previous case, we return a 1d vector result rather than a 2d array with 1 row and ncolumns.
Finally we have the case of
\[\left[\begin{array}{cc} 0 & 1\ 0 & 0\end{array}\right] \left[\begin{array}{cc} 0 & 0\ 1 & 0\end{array}\right] = \left[\begin{array}{cc} 1 & 0\ 0 & 0\end{array}\right] \]
(matrix-multiply [[0 1]
[0 0]]
[[0 0]
[1 0]])✓
This example is adapted from here. The correct answer is at the bottom of that page, and shown here.
\[\left[\begin{array}{cccc} 1 & 2 & -2 & 0 \ -3 & 4 & 7 & 2 \ 6 & 0 & 3 & 1\end{array}\right] \left[\begin{array}{cc} -1 & 3 \ 0 & 9 \ 1 & -11 \ 4 & -5 \end{array}\right] = \left[\begin{array}{cc} -3 & 43 \ 18 & -60 \ 4 & -5\end{array}\right] \]
For readability we use temporary variables here, and pretty-print the result.
(let ((A [[1 2 -2 0]
[-3 4 7 2]
[6 0 3 1]])
(B [[-1 3]
[0 9]
[1 -11]
[4 -5]]))
(pp (matrix-multiply A B)))✓
So, all these example work as we expect. The elisp function for matrix-multiply does a lot of work for you to make these cases work, including error checking for dimensional consistency.
It was not any easier to write this dynamic module than the previous one I used with the Gnu Scientific Library. The approach and code are remarkably similar. In one way the GSL was easier to use; it worked out of the box, whereas I had to fiddle with a security option in my OS to get it to run MKL! My Anaconda Python distribution must get around that somehow since it ships with an MKL compiled Numpy and scipy.
The idea of using elisp for analysis of the inputs and making sure they are correct is a good one and helps prevent segfaults. Of course it is a good idea to write defensive c-code to avoid that too. Overall, this is another good example of expanding the capabilities of Emacs with a dynamic module.
The c-code is loosely adapted from https://software.intel.com/en-us/node/529735. We do not implement the full dgemm behavior which is able to calculate
#include "emacs-module.h"
#include "emacs-module-helpers.h"
#include <mkl.h>
int plugin_is_GPL_compatible;
emacs_value Fmkl_dgemm (emacs_env *env, ptrdiff_t nargs, emacs_value args[], void *data)
{
double *A, *B, *C;
int m, n, k, i, j;
double alpha = 1.0;
double beta = 0.0;
// These will be 2d vectors
emacs_value M0 = args[0]; // array 1 - A (m x k)
emacs_value M1 = args[1]; // array 2 - B (k x n)
// I need to get the number of rows and columns of each one.
m = env->vec_size(env, M0);
k = 0;
// We assume that we have a 2d array.
emacs_value el1 = env->vec_get (env, M0, 0);
k = env->vec_size(env, el1);
// Now we know A has dimensions (m, k)
emacs_value el2 = env->vec_get (env, M1, 0);
n = env->vec_size(env, el2);
// Now we know M1 had dimensions (k, n)
// Now we have to build up arrays.
// We are looking at a * b = c
A = (double *)mkl_malloc( m*k*sizeof( double ), 64 );
B = (double *)mkl_malloc( k*n*sizeof( double ), 64 );
C = (double *)mkl_malloc( m*n*sizeof( double ), 64 );
if (A == NULL || B == NULL || C == NULL) {
printf( "\n ERROR: Can't allocate memory for matrices. Aborting... \n\n");
mkl_free(A);
mkl_free(B);
mkl_free(C);
return 1;
}
//populate A
emacs_value row, ij;
for (int i = 0; i < m; i++)
{
row = env->vec_get(env, M0, i);
for (int j = 0; j < k; j++)
{
// get M0[i, j]
ij = env->vec_get(env, row, j);
A[k * i + j] = extract_double(env, ij);
}
}
// populate B
for (int i = 0; i < k; i++)
{
row = env->vec_get(env, M1, i);
for (int j = 0; j < n; j++)
{
// get M0[i, j]
ij = env->vec_get(env, row, j);
B[n * i + j] = extract_double(env, ij);
}
}
// initialize C. The solution will have dimensions of (rows1, cols2).
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
C[n * i + j] = 0.0;
}
}
// the multiplication is done here.
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
m, n, k, alpha, A, k, B, n, beta, C, n);
// now we build up the vector to return
emacs_value vector = env->intern(env, "vector");
emacs_value *array = malloc(sizeof(emacs_value) * m);
emacs_value *row1;
emacs_value vec;
for (int i = 0; i < m; i++)
{
row1 = malloc(sizeof(emacs_value) * n);
for (int j = 0; j < n; j++)
{
row1[j] = env->make_float(env, C[j + i * n]);
}
vec = env->funcall(env, vector, n, row1);
array[i] = vec;
free(row1);
}
emacs_value result = env->funcall(env, vector, m, array);
free(array);
return result;
}
int emacs_module_init(struct emacs_runtime *ert)
{
emacs_env *env = ert->get_environment(ert);
DEFUN("mkl-dgemm", Fmkl_dgemm, 2, 2,
"(mkl-dgemm A B)\n"\
"Multiply the matrices A and B. A and B must both be 2d vectors.\n" \
"Returns the product as a vector.",
NULL);
provide(env, "mkl");
return 0;
}To build this we have to run elisp:org-babel-tangle to generate the mkl.c file, and then run this shell block to compile it.
sh /opt/intel/mkl/bin/mklvars.sh intel64
gcc -Wall -m64 -I${MKLROOT}/include -fPIC -c mkl.c
gcc -shared -L${MKLROOT}/lib -Wl,-rpath,${MKLROOT}/lib -lmkl_rt -lpthread -lm -ldl -L. -lemacs-module-helpers -o mkl.so mkl.oWe will often want to know the shape of our arrays. The shape is how many elements there are in each dimension. Here we define a recursive function that gets the shape of arbitrarily nested vectors and returns a vector of the shape. We define some helper functions to get the number of dimensions, elements, rows and columns.
The main function is a helper elisp function that multiplies two arrays. The function analyzes the shapes and transforms 1d vectors into the right 2d shape to multiply them together, and then returns the shape that makes sense. The c-code is not very robust to mistakes in the array dimensions. It tends to make emacs segfault if you get it wrong. So we try to avoid that if possible.
We have four cases to consider for multiplication:
- 2d * 2d
- (assert (= m1 n2)) return (n1, m2)
- 1d * 2d
- 1d is a row vector (1, n1) (assert (= n1 m2)) return vector with n2 elements
- 2d * 1d
- 1d is a column vector (m2, 1) (assert (= n1 m2)) return vector with m2 elements
- 1d * 1d
- (assert (= (length m1) (length m2)) return a scalar
Here is the
(add-to-list 'load-path (expand-file-name "."))
(require 'mkl)
(require 'cl)
(require 'seq)
(defun vector-shape (vec)
"Return a vector of the shape of VEC."
(let ((shape (vector (length vec))))
(if (vectorp (aref vec 0))
(vconcat shape (vector-shape (aref vec 0)))
shape)))
(defun vector-ndims (vec)
"Returns the number of dimensions in VEC."
(length (vector-shape vec)))
(defun vector-numel (vec)
"Returns the number of elements in VEC."
(if (> (length vec) 0)
(seq-reduce '* (vector-shape vec) 1)
0))
(defun vector-nrows (vec)
"Return the number of rows in VEC."
(aref (vector-shape vec) 0))
(defun vector-ncols (vec)
"Return the number of columns in VEC."
(aref (vector-shape vec) 1))
(defun matrix-multiply (A B)
"Return A * B in the matrix multiply sense."
(cond
;; 1d * 1d i.e. a dot-product
((and (= 1 (vector-ndims A))
(= 1 (vector-ndims B))
(= (length A) (length B)))
;; this is easy to compute so we don't use dgemm.
(seq-reduce '+ (mapcar* (lambda (a b) (* a b)) A B) 0.0))
;; 2d * 1d (m1, n1) * (n2, 1)
((and (= 2 (vector-ndims A))
(= 1 (vector-ndims B))
;; ncols-A = len-B
(= (vector-ncols A) (length B)))
;; transform B into a 2d column vector
(let* ((B2d (apply 'vector (mapcar 'vector B)))
(result (mkl-dgemm A B2d)))
;; Now call (dgemm A B2d) -> (m2, 1) column vector
;; and convert it back to a 1d result
(cl-map 'vector (lambda (v) (aref v 0)) result)))
;; 1d * 2d (1, n1) * (m2, n2) len-A = nrows-B
((and (= 1 (vector-ndims A))
(= 2 (vector-ndims B))
(= (length A) (vector-nrows B)))
;; transform B into a 2d row vector
(let* ((A2d (vector A))
(result (mkl-dgemm A2d B)))
;; should be a 2d row vector
(aref result 0)))
;; 2d * 2d (m1, n1) * (m2, n2) rows-A = ncols-B
((and (= 2 (vector-ndims A))
(= 2 (vector-ndims B))
(= (vector-ncols A)
(vector-nrows B)))
;; call (dgemm A B) and return result
(mkl-dgemm A B))
(t
;; Error checking, getting here means none of the cases above were caught.
;; something is probably wrong.
(cond
((or (> (vector-ndims A) 2)
(> (vector-ndims B) 2))
(error "One of your arrays has more than 2 dimensions. Only 1 or 2d arrays are supported"))
((and (= 1 (vector-ndims A))
(= 1 (vector-ndims B))
(not (= (length A) (length B))))
(error "A and B must be the same length.
len(A) = %d
len(B) = %d" (length A) (length B)))
((and
(= (vector-ndims A) 2)
(= (vector-ndims B) 2)
(not (= (vector-nrows A) (vector-ncols B))))
(error "Your array shapes are not correct.
The number of rows in array A must equal the number of columns in array B.
There are %d rows in A and %d columns in B" (vector-nrows A) (vector-ncols B)))
((and
(= (vector-ndims A) 2)
(= (vector-ndims B) 1)
(not (= (vector-nrows A) (length B))))
(error "Your array shapes are not correct.
The number of rows in array A must equal the number of columns in array B.
There are %d rows in A and %d columns in B" (vector-nrows A) (length B)))
(t
(error "Unknown error"))))))