7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ If your matrix is sparse, you might consider using MKL's PARDISO, which will provide both left and right solutions (cf. examples in the directory /opt/intel/compilers_and_libraries_2017/<YOUR_PLATFORM_HERE>/mkl/examples/solverc/source/). $\endgroup$ – NewDogOldTricks Nov 1 '17 at 22:08
9
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ That is something I have tried before, but for some reason this method (implemented in octave) returns an empty array for a sparse (mostly diagonal) matrix. Is there anything you can possibly suggest about this particular C implementation? $\endgroup$ – ivan-k Dec 2 '13 at 1:07
  • $\begingroup$ It's not clear what you tried in octave, since as I described it, it should never return an empty matrix. You simply cannot use xGESV on non-square matrices, period. $\endgroup$ – Victor Liu Dec 2 '13 at 1:15
  • $\begingroup$ Fair enough. Apologies for unclear question. I'l tune back in after I have tried it. Thanks again. $\endgroup$ – ivan-k Dec 2 '13 at 1:29
  • $\begingroup$ Both approaches work well for smaller matrices that I used for testing. I know that I might be asking too much, but kindly consider my last question. Can you comment on numerical stability and the speed for the two approaches (QR vs SVD)? Thank you. $\endgroup$ – ivan-k Dec 2 '13 at 22:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.