topic How to use mkl to solve large scale sparse least square prb? in Intel® oneAPI Math Kernel Library

How to use mkl to solve large scale sparse least square prb?

liunsteingmail_com — Fri, 19 Jun 2009 11:51:31 GMT

Hi all:

I understand this is more of a math problem than a software one
But i do hope some people can give me some hints

My problem is of the type classic: Ax=b and A is my Jacobian matrix(the orginal least square prb is nonlinear, i use gauss-newton method)
and A is sparse and rectangle
for the moment i use the normal equations: A'Ax=A'b
and for this i use functions in sparse blas
and then i solve the resulting system with pardiso

i'm sure there are some better ways to do this
for the moment i'm still struggling to make this most 'simple' method work

thanks in advance

Re: How to use mkl to solve large scale sparse least square prb

ArturGuzik — Sat, 20 Jun 2009 06:50:03 GMT

Quoting - liunsteingmail.com

My problem is of the type classic: Ax=b and A is my Jacobian matrix(the orginal least square prb is nonlinear, i use gauss-newton method)
and A is sparse and rectangle
for the moment i use the normal equations: A'Ax=A'b
and for this i use functions in sparse blas
and then i solve the resulting system with pardiso

i'm sure there are some better ways to do this
for the moment i'm still struggling to make this most 'simple' method work

thanks in advance

Hi,

if youre problem is non-linear and you try to use normal equations approach in between, it sounds like asking for trouble. The're two approaches for solving linear LS problem, one based on QR decomposition and the other using SVD (for more tricky things). See Numerical Recipes for a discussion on that (I believe it's in Chapter on data modelling). For a start you can just try MKL routines for that (LAPACK for LS) which offer (copying from the docs):

bd bidiagonal matrix
ge general matrix
gb general band matrix
hs upper Hessenberg matrix
or (real) orthogonal matrix
op (real) orthogonal matrix (packed storage)
un (complex) unitary matrix
up (complex) unitary matrix (packed storage)
pt symmetric or Hermitian positive-definite tridiagonal matrix
sy symmetric matrix
sp symmetric matrix (packed storage)
sb (real) symmetric band matrix
st (real) symmetric tridiagonal matrix
he Hermitian matrix
hp Hermitian matrix (packed storage)
hb (complex) Hermitian band matrix
tr triangular or quasi-triangular matrix.

But, what you actually should try is MKL reverse communication, to try to manage your original (non-linear) problem. Is the sparsity that significant it's worth going through this?

A.

Re: How to use mkl to solve large scale sparse least square prb

liunsteingmail_com — Sun, 21 Jun 2009 16:46:23 GMT

Hi thanks for the reply
I didn't consider lapack for a LS solution since the matrix is really huge and i have to take advantage of the sparsity. otherwise the problem will be intracable.

I have read some theories on LS problem. espetially the classic reference books from SIAM. but it's another story putting these theories in pratice. actually i would just like an outline of the solution. things like what are the major steps when dealing with this problem. and what are the mkl functions people use to solve them.

I have great difficulties understanding the parameters in mkl functions even with the lengthy user manual. the meaning of most of the function arguments i can not find in matrix computation books. so in the end i don't know how to choose them.

Re: How to use mkl to solve large scale sparse least square prb

ArturGuzik — Mon, 22 Jun 2009 00:04:34 GMT

I have great difficulties understanding the parameters in mkl functions even with the lengthy user manual. the meaning of most of the function arguments i can not find in matrix computation books. so in the end i don't know how to choose them.

With parameters should no problem. MKL people can always respond to any related issues. That's for sure. I believe your main concern should be the non-linearity of the problem.

A.

Re: How to use mkl to solve large scale sparse least square prb

ArturGuzik — Mon, 22 Jun 2009 00:15:26 GMT

I had once a book on that. It listed software packages and some implementation details. See here.

Dis you try these routines?

From the comments inside:

!-------------------------------------------------------------------
! lstp generates a sparse least-squares test problem of the form
! ( A )*x = ( b )
! ( damp*I ) ( 0 )
! for solution by LSQR, or a sparse underdetermined system
! Ax + damp*s = b
! for solution by LSQR or CRAIG. The matrix A is m by n and is
! constructed in the form A = Y*D*Z, where D is an m by n
! diagonal matrix, and Y and Z are Householder transformations.
!
! m and n may contain any positive values.
! If m >= n or damp = 0, the true solution is x as given.
! Otherwise, x is modified to contain the true solution.

Anyway, If you decide to stuck with MKL, and you have some difficulty in calling MKL I'm sure somebody can assist and make code working. You just need to post more details.

A.