general linear model example -- constrained least squares using foreign function calls to the LAPACK library

In this example, we make foreign function calls to LAPACK to solve constrained least squares problems. LAPACK is already used by Macaulay2 for a number of its linear algebra computations, but some functions aren't available.

One of these is dggglm, which solves constrained least squares problems that have applications to general Gauss-Markov linear models in statistics. In particular, suppose $A$ is an $n\times m$ real matrix, $B$ is an $n\times p$ real matrix, and $\mathbf d$ is a vector in $\RR^n$. We would like to find the vectors $\mathbf x\in\RR^m$ and $\mathbf y\in\RR^p$ that minimize $\|\mathbf y\|$ subject to the constraint $\mathbf d = A\mathbf x + B\mathbf y$.

The dggglm function does exactly this. Its signature is void dggglm_(int *n, int *m, int *p, double *A, int *ldA, double *B, int *ldB, double *d, double *x, double *y, double *work, int *lwork, int *info). Note the trailing underscore that was added by the Fortran compiler. Each of the 13 arguments are pointers, so we use voidstar to represent them in our call to foreignFunction. Since Macaulay2 is already linked against LAPACK, we don't need to worry about calling openSharedLibrary.

The parameters n, m and p are exactly the numbers $n$, $m$, and $p$ above. The parameters A, B, and d are arrays that will store the entries of $A$, $B$, and $\mathbf d$. LAPACK expects these to be in column-major order, in contrast to the row-major order used by Macaulay2. So we add helper methods to take care of this conversion. The local variable T is a ForeignArrayType created by ZZ * ForeignType.

The parameters ldA and ldB are the "leading dimensions" of the arrays A and B. We will use $n$ for both. The arrays x and y will store the solutions. The array work has dimension lwork (for which we will use $n + m + p$) that the procedure will use as a workspace. Finally, info stores the return value (0 on success).

The following method prepares matrices $A$ and $B$ and a vector $\mathbf d$ to be sent to the foreign function we created above, and then converts the solutions into a sequence two of Macaulay2 vectors. Note that we use getMemory to allocate memory for x, y, and info (which we name i here to avoid conflicting with info) and address to get pointers to the other integer arguments. We also use ForeignType * voidstar to dereference i and determine whether the call was successful.

Finally, let's call this method and solve a constrained least squares problem. This example is from Kourouklis, S., & Paige, "A Constrained Least Squares Approach to the General Gauss-Markov Linear Model".