Matrix multiplication and inverse problems with accelerate framework

Question

I am trying to multiply two matrices in Objective-C. I have imported the accelerate framework to my project in Xcode, everything compiles just fine. I did the matrix multiplication on my calculator and got correct values, however, when running the code, I do not. This is how I defined my matrices..

float matrixA [3][3] = {{-.045, -.141, -.079}, {-.012, -.079, .0578}, {.112, -.011, -.0830}};

float matrixB [3][1] = {{40}, {-61}, {98}};

I then used the mmul function in the accelerate framework..

vDSP_mmul(&matrixA[3][3], 1, &matrixB[3][1], 1, &results[3][1], 1, 3, 1, 3);

The array results was created by doing the following..

float results[3][1];

I just placed this in the viewDidLoad method of an empty project so I could NSLog the results. So when I multiply matrixA by matrixB I should get the following results.. (-1, 10, -3). However, the results in the NSLog show (-0.045, 0.000, 0.000). This is not correct and I don't understand why. My understanding was that this function would multiply two matrices together. But I am not sure what it is doing. I am probably inputing something incorrectly and was hoping someone could help me out.

Side Note: matrixA is really the inverse of another matrix. However, I can't find anything in the accelerate framework that will take the inverse. I did find a function called

sgetrf_

with LAPACK but don't really get it. If anyone has any help, advice, or some tutorial to follow I would appreciate it, been at this for three days now looking thing up on the internet!!

Answer 1

Let's walk through three different ways to do this computation (including the inverse) on OS X or iOS.

First, let's do what you were more-or-less trying to do to start with: we'll do the computation using LAPACK and BLAS from the Accelerate framework. Note that I've used the BLAS function cblas_sgemv instead of vDSP_mmul to do the matrix-vector multiplication. I've done this for three reasons. First, it's more idiomatic in code that uses LAPACK. Second, LAPACK really wants matrices stored in column-major order, which BLAS supports, but vDSP does not. Finally, vDSP_mmul is actually just a wrapper around the BLAS matrix-multiply routines, so we can cut out the middleman. This will work on OS X back to 10.2 and iOS back to 4.0--i.e. on any target you can reasonably expect to encounter today.

#include <Accelerate/Accelerate.h>
#include <stdio.h>
#include <string.h>

void using_blas_and_lapack(void) {
    printf("Using BLAS and LAPACK:\n");
    //  Note: you'll want to store A in *column-major* order to use it with
    //  LAPACK (even though it's not strictly necessary for this simple example,
    //  if you try to do something more complex you'll need it).
    float A[3][3] = {{-4,-3,-5}, {6,-7, 9}, {8,-2,-1}};
    float x[3] = { -1, 10, -3};
    //  Compute b = Ax using cblas_sgemv.
    float b[3];
    cblas_sgemv(CblasColMajor, CblasNoTrans, 3, 3, 1.f, &A[0][0], 3, x, 1, 0.f, b, 1);
    printf("b := A x = [ %g, %g, %g ]\n", b[0], b[1], b[2]);
    //  You probably don't actually want to compute A^-1; instead you simply
    //  want to solve the equation Ay = b for y (like y = A\b in matlab).  To
    //  do this with LAPACK, you typically use the routines sgetrf and sgetrs.
    //  sgetrf will overwrite its input matrix with the LU factorization, so
    //  we'll make a copy first in case you need to use A again.
    float factored[3][3];
    memcpy(factored, A, sizeof factored);
    //  Now that we have our copy, go ahead and factor it using sgetrf.
    __CLPK_integer n = 3;
    __CLPK_integer info = 0;
    __CLPK_integer ipiv[3];
    sgetrf_(&n, &n, &factored[0][0], &n, ipiv, &info);
    if (info != 0) { printf("Something went wrong factoring A\n"); return; }
    //  Finally, use the factored matrix to solve y = A\b via sgetrs.  Just
    //  like sgetrf overwrites the matrix with its factored form, sgetrs
    //  overwrites the input vector with the solution, so we'll make a copy
    //  before calling the function.
    float y[3];
    memcpy(y, b, sizeof y);
    __CLPK_integer nrhs = 1;
    sgetrs_("No transpose", &n, &nrhs, &factored[0][0], &n, ipiv, y, &n, &info);
    printf("y := A\\b = [ %g, %g, %g ]\n\n", y[0], y[1], y[2]);
}

Next, let's do the computation using LinearAlgebra.h, which is a new feature of the Accelerate Framework in iOS 8.0 and OS X 10.10. It abstracts away nearly all of the bookkeeping necessary for these computations, but it only offers a tiny piece of the full functionality that available in LAPACK and BLAS. Note that while there's some boilerplate necessary to convert our raw arrays to and from la_objects, once we have them the actual computation is very straightforward.

#include <Accelerate/Accelerate.h>
#include <stdio.h>

void using_la(void) {
    printf("Using LA:\n");
    //  LA accepts row-major as well as column-major data, but requires a short
    //  two-step dance to get our data in.
    float Adata[3][3] = {{-4, 6, 8}, {-3,-7,-2}, {-5, 9,-1}};
    float xdata[3] = { -1, 10, -3};
    la_object_t A = la_matrix_from_float_buffer(&Adata[0][0], 3, 3, 3, LA_NO_HINT, LA_DEFAULT_ATTRIBUTES);
    la_object_t x = la_vector_from_float_buffer(xdata, 3, 1, LA_DEFAULT_ATTRIBUTES);
    //  Once our data is stored as LA objects, it's easy to do the computation:
    la_object_t b = la_matrix_product(A, x);
    la_object_t y = la_solve(A, b);
    //  And finally we need to get our data back out:
    float bdata[3];
    if (la_vector_to_float_buffer(bdata, 1, b) != LA_SUCCESS) {
        printf("Something went wrong computing b.\n");
        return;
    } else printf("b := A x = [ %g, %g, %g ]\n", bdata[0], bdata[1], bdata[2]);
    float ydata[3];
    if (la_vector_to_float_buffer(ydata, 1, y) != LA_SUCCESS) {
        printf("Something went wrong computing y.\n");
        return;
    } else printf("y := A\\b = [ %g, %g, %g ]\n\n", ydata[0], ydata[1], ydata[2]);
}

And finally, one more method. If your matrices are really just 3x3, this is the method I would use, because the overhead involved in either of the previous methods will swamp the actual computation. However, this only works for matrices of size 4x4 or smaller. There's another new header in iOS 8.0 and OS X 10.10 specifically focused on small matrix and vector math, that makes this really simple and efficient:

#include <simd/simd.h>
#include <stdio.h>

void using_simd(void) {
    printf("Using <simd/simd.h>:\n");
    matrix_float3x3 A = matrix_from_rows((vector_float3){-4,  6,  8},
                                         (vector_float3){-3, -7, -2},
                                         (vector_float3){-5,  9, -1});
    vector_float3 x = { -1, 10, -3 };
    vector_float3 b = matrix_multiply(A, x);
    printf("b := A x = [ %g, %g, %g ]\n", b[0], b[1], b[2]);
    vector_float3 y = matrix_multiply(matrix_invert(A),b);
    printf("y := A^-1 b = [ %g, %g, %g ]\n\n", y[0], y[1], y[2]);
}

And finally, let's just double-check that these all give identical results (up to small rounding differences):

scanon$ xcrun -sdk macosx clang matrix.m -framework Accelerate  && ./a.out
Using BLAS and LAPACK:
b := A x = [ 40, -61, 98 ]
y := A\b = [ -0.999999, 10, -3 ]

Using LA:
b := A x = [ 40, -61, 98 ]
y := A\b = [ -1, 10, -3 ]

Using <simd/simd.h>:
b := A x = [ 40, -61, 98 ]
y := A^-1 b = [ -1, 10, -3 ]

Answer 2

You are passing pointers to memory after the end of your matrices.

Fixing that would look like this (untested code):

vDSP_mmul(&matrixA[0][0], 1, &matrixB[0][0], 1, &results[0][0], 1, 3, 1, 3);

Passing an array to a function in C will actually pass a pointer to the first element of the array, which in this case appears to be what you need to do. You were passing a pointer to a piece of memory directly after the final array element in your matrices, which means you'll get nonsense results, or a crash.