Solving for null space of a matrix with mkl LAPACK

Question

I want to find a solution for $xA=0$, where $A$ is a square matrix. I know that most of the LAPACK routines solve for $Ax=b$. So I take $A^T$ as a, and set $b=0$. I have an additional restriction of $\Sigma x=1$.Thus, I add an additional row of $ones$ to $A^T$ and a single $one$ to $b$. This also allows to avoid the trivial solution in which $x=0$.

This approach works fine in Octave, but I am struggling to implement this with LAPACK. I keep getting info=4, which means that my fourth factor $U$ in the $LU$ factorization of $A^T$ with a row of $ones$ is singular. This is happening because my last row is $ones$.

Is there a good way to solve $xA=0$? I should add that in my real data, $A$ is sparse.

Here is my failing attempt in C, where b=[0,0,0,1].

float *A, *b;
/* Ax = b
 * Ax = xT(A)   //transpose
 * A is m x m   //matrix   
 * b is m x 1  //vector
 * x is m x 1   //vector
 */
int m = 3;
int scale = 1;
A = (float *)mkl_malloc(m*m,32);
A[0]=-5;A[1]= 2;A[2]= 3;
A[3]= 4;A[4]= -10;A[5]= 6;
A[6]= 7;A[7]=8;A[8]= -15;

b = (float*)mkl_malloc(m+1,32);
for (int i = 0; i < m; i++) {
    b[i]=0.0;
}
b[m]=1.0;

int matrix_order = LAPACK_ROW_MAJOR;

int ipiv[m+1];

int info;
//transpose and add a row of ones
float * Acopy = (float*)mkl_malloc((m+1)*m,32);
mkl_somatcopy('R','T',m,m,scale,A,m,Acopy,m);
//assign last row of ones
for (int i = 0; i < m; i++){
    int p = m*m+i;
    Acopy[p]=1.0;
}

printf("Original\n");
print_matrix(A, m, m );
printf("Transpose with ones\n");
print_matrix(Acopy, m+1, m );

info = LAPACKE_sgesv(matrix_order,m+1,1,Acopy,m+1,ipiv,b,1);

printf("$?=%d\n",info);

print_matrix(b, m+1, 1 );
return 0;

Thank you.

Answer 1

There's no reason to append a row of 1's. You should just perform a rank-revealing QR factorization (like with routine SGEQP3) on $A^T$, and the last column of $Q$ should be in the nullspace. This has the added advantage that the relative magnitude of the last element on the diagonal of $R$ gives you some idea of how singular the solution is. Even better would be to perform a singular value decomposition of $A=U\Sigma V^T$. Then, the last column of $U$ is a singular solution to $xA=0$. You can simply normalize to achieve $\sum x = 1$ afterwards.

Note that you can't just append a row of ones at the bottom of $A$ since then you have a non-square system. You can probably solve that by least squares using SGELS (or any of the other variants). This generally is less robust than the previous two suggestions I gave, although it enforces the normalization explicitly.

Edit: Note on speed and stability. I would guess QR is slightly faster than SVD, while SVD is more numerically stable (although practically, they are sufficiently stable that people don't worry about it). I don't have a good sense of speed since more sophisticated pivoting in QR is slower. Note finally that these approaches only work for dense matrices. If you want something that works in the sparse case, you can try power iteration on the inverse of $A$. This should work as long as $A$ is not exactly singular (which basically never happens numerically).