3

Solve Ax = b. Real double. A is overdetermined Mx2 with M >> 2. b is Mx1. I've run a ton of data against mldivide, and the results are excellent. I wrote a mex routine with MKL LAPACKE_dgels and it's nowhere near as good. The results have a ton of noise and the underlying signal is barely there. I checked the routine against the MKL example results first. I've searched through the mldivide doc (flowchart) and the SO questions. All I found is Matlab uses QR factorization for overdetermined rectangular.

What should I try next? Am I using the wrong LAPACK routine? Please help guide me in the right direction.

Update: To within E-15 floating point difference on the solution vector, Intel MKL LAPACKE_dgels has the same result as Matlab mldivide for real double overdetermined (rectangular) problems. As far as I can tell, this is the QR method used.

Beware the residuals returned from this dgels. They do not equate to b - Ax. Many of them are close to this value, whereas some are far from it.

4
  • For us to check if you are using the function the wrong way, you'd better provide a minimum working example (even if with small matrices the issue is not reproduced).
    – Mefitico
    Commented Sep 30, 2019 at 16:11
  • As an alternative, give a try to LAPACK's DGELSY(). It also makes use of QR decomposition to solve overdetermined systems. Finally, give a try to DGELSD() if the matrix A is ill conditionned. See stackoverflow.com/questions/41637108/… and stackoverflow.com/questions/55367024/…
    – francis
    Commented Sep 30, 2019 at 18:43
  • 1
    @francis. Thank you, that's good advice. I read Armadillo uses QR then falls back to SVD (dgelsd) if its ill conditioned. But the int info return parameter of dgels is never positive (always zero). Commented Sep 30, 2019 at 18:56
  • Indeed, the file include/armadillo_bits/glue_solve_meat.hpp describes a call to status = auxlib::solve_approx_fast() and a call to auxlib::solve_approx_svd() if status is false. solve_approx_fast() calls lapack::gels and auxlib::solve_approx_svd() calls lapack::gelsd(). The return parameter of dgels become positive if a null pivot is found, thus preventing the system from being solved. The fact that it is always zero is consistent with a good conditioning of your matrix. Congratulations for having solved your problem yourself!
    – francis
    Commented Sep 30, 2019 at 20:56

1 Answer 1

5

The problem was not the solution x, rather the returned residuals from DGELS. This routine's outputs are modify-in-place on the input array pointers. The MKL doc says the input array b is overwritten with the output vector x for the first N rows, then the residuals in N+1 to M. I confirmed this with my code.

The mistake was in aligning the b[N+1] residuals to original inputs b[1], and making further algorithmic decisions on that. The correct alignment of residual to original input is b[1] to b[1]. The first N residuals are not available; you have to compute those afterwards.

The doc doesn't say they are residuals per se, rather specifically

the residual sum of squares for the solution in each column is given by the sum of squares of modulus of elements n+1 to m in that column.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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