Cholesky factorization, traspose and inversion of sparse matrix.

[Fiori]

Hello!

Let Sigma be a sparse matrix. I would like to compute the Cholesky factorization of Sigma (the Upper(Lt) or lower triangular (L)), transpose it, and compute the folowing terms

w = inv(L)*mu;
m = inv(Lt)*w;
v = inv(Lt)*b;

where mu, b are known.

The problem I face is that I can't find the routines (and examples) when the matrix is sparse. Any help would be really appreciated.
Thank you very much.

[Ying_H_Intel]

Hello, 

What are your problem size and sparse matrix structure (like Non-zero number ratio, storage format) etc? 

I guess, the sparse solver and sparse BLAS in MKL may help.  Please refer to the manual of MKL  or

 https://software.intel.com/en-us/articles/intel-mkl-sparse-solvers-training-material  ( i will attach the pdf file there). 

and you can find code sample in MKL install forder : solverc or solverf (fortran). 

Sparse Solver Routines=> Intel MKL PARDISO

Best Regards,

Ying 

[mecej4]

If all that is required is to solve A.x = b given a sparse matrix A and vectors x and b, it is neither necessary nor efficient to form and use the explicit inverse of A. When you use a sparse linear solver, the usual sequence of operations is (i) analyze (ii) factorize and (iii) solve. If the matrix is banded, you can use the band solver routines in MKL/Lapack. Otherwise, use the sparse solver Pardiso directly or through the DSS interface. For details, consult the MKL documentation and the example source and data files in the MKL examples/solverf directory.

[Fiori]

Hello. Thank you very much for your help.

I want to sample from a Multivariate Normal distribution(Q^{-1}b, Q^{-1}) , where Q is a triagonal, symmetric matrix:

the elements of the main diagonal are {a1, a2, a2, ...., a2, a1} and the elements of the upper (down) diagonal are {a3, a3, a3, .... a3}.

The algorithm in order to do this is the following:

1. Compute the Cholesky factorization, Q=LU, where U=L^{T}

2. Solve Lw = b

3. Solve Um=w

4. Sample z~N(0,1)

5. Solve Uv=z

6.Compute x=m+y

7.Return x

The number of rows of Q is 5e+5. 

Steps 2 and 3 gives the solution of Qm=b. The function PARDISO computes it. How can I compute step 5? Is it better to write it as a band or as a sparse matrix?   

[Ying_H_Intel · Accepted solution]

Hi GF,

As I understand, both Lapack function with band format and pardiso can do your computation. You may try both.

But if your matrix format is band obviously, I evaluate a  band matrix and lapack function call may better (easy to use and better performance). 

 Pardiso support separate Forward and Backward Substitution, so you should be able to call phase 331 for different & multiply z . 

The solver execution step (see parameterphase= 33below) can be divided into two or three separate

substitutions: forward, backward, and possible diagonal. This separation can be explained by the examples of
solving systems with different matrix types.
A real symmetric positive definite matrix A(mtype= 2) is factored by Intel MKL PARDISO as A= L*L
T
. In
this case the solution of the system A*x=bcan be found as sequence of substitutions: L*y=b(forward
substitution, phase=331) andL
T
*x=y(backward substitution, phase=333)

If iparm[35] = 2, phases 331, 332, and 333 perform a different

decomposition:
You can supply a custom implementation for phase 332 instead of calling
pardiso. For example, it can be implemented with dense LAPACK
functionality. Custom implementation also allows you to substitute the
matrix Swith your own.
NOTE
For very large Schur complement matrices use LAPACK functionality to
compute the Schur complement vector instead of the Intel MKL
PARDISO phase 332 implementation.

Best Regards,

Ying