lapack

HI Pavel,

               I believe MAtrix Toolkit Java already uses Lapack for computing Cholesky. Unfortunately, it still takes quite some time to compute.

Nothing is faster than lapack in terms of algorithms. It's just Java wrappers that are slow.

Try to solve D=A^-1*M*(A^T)^-1 across iterations

@ Konstantin , I realized that I was missing the unknown variable 'X'  in my original post.  Also, do note that A is a rectangular matrix

Jayanth^

try some packages from www.netlib.org searching word "sparse"

http://www.netlib.org/cgi-bin/search.pl

for example:

http://www.netlib.org/sparse/

A is rectangula - so you are to minimize ||A*D*A^T-M||^2. May be with the help of iterative Broyden–Fletcher–Goldfarb–Shanno algorithm

https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm

Jayanth,

First I would like to highlight the answer of Fasino. That will really work out!

For the sake of curiosity, let me suggest an alternative solution.

Let v_i denote the i-th column of, that is A = [v1 v2 ... vn]. At the (k+1)-th iteration,

you already have the solution x^{(k)} of PHI_k x = b,where PHI_k = A D_k A^T.

Denote by x^{(k+1)} the solution of of PHI_k+1 x = b,where

PHI_k+1 = A D_k+1 A^T. (1)

It is interesting to consider the new diagonal matrix D_k+1 as a perturbation of D_k, that is D_k+1 = D_k + delta.

Let us denote by d_i, i = 1, 2, ..., m, the diagonal elements of delta.

Then, we can rewrite the PHI_k+1 in (1) as

PHI_k+1 = PHI_k + sum_{i=1}^n d_i v_i v_i^T.

In order to obtain x^{(k+1)}, let us introduce a new sequence of vectors defined by

z^{(0)} := x^{(k)} and

z^{(j)} is the solution of [ PHI_k + sum_{i=1}^j d_i v_i v_i^T ] z = b, j 1, 2, ..., n.

Therefore, Our target is to compute z^{(n)} = x^{(k+1)}.

Set B_0 = PHI_k, B_j = PHI_k + sum_{i=1}^j d_i v_i v_i^T, j = 1, 2, ..., n, so that

B_j z^{(j)} = b , j = 0, 1, 2, ..., n.  (2)

The important thing here is that, to pass from B_j to B_j+1, we add the 1-rank matrix d_j+1 v_j+1 v_j+1^T.

Therefore, the new solution z^{(k+1)} can be easily obtained from z^{(k)} in terms of the Sherman-Morisson formula

https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula

More precisely, we have

z^{(j+1)} = [B_{j+1}]^-1 b = [B_{j} + d_j+1 v_j+1 v_j+1^T]^-1 b

= [B_j - \frac{d_j+1}{1 + d_j+1 *(v_j+1^T B_j^-1 v_j+1)} B_j^-1 v_j+1 v_j+1^T B_j^-1 ] b

= z^{(j)} - \frac{d_j+1}{1 + d_j+1 *(v_j+1^T B_j^-1 v_j+1)}[B_j^-1 v_j+1 v_j+1^T B_j^-1] b

= z^{(j)} - \frac{d_j+1}{1 + d_j+1 *(v_j+1^T B_j^-1 v_j+1)}[B_j^-1 v_j+1 v_j+1^T] z^{(k)} (3)

The laborious part in the right-hand side of (3) is the calculus of B_j^-1 v_j+1. Nevertheless, notice that

B_j^-1 v_j+1 is the solution of B_j x = v_j+1, which is exactly the system (2), but with a different independent

term, namely b' = v_j+1. Hence, we can use an iteration process as above in order to precompute these terms.

For instance, in the computation of B_1^-1 b', we will need to compute B_0^-1 b' in (3).

In the computation of B_2^-1 b', we use the term B_1^-1 b' which we already computed and so on.

Remark1: The vectors B_0^-1 v_1, B_0^-1 v_2, ... , B_0^-1 v_n can be computed only once and then stored (notice that they are independent of k). For each value of j, only one vector B_j^-1 v_j+1 is required and it can be computed on the fly from the previous one at  cost of O(m) operations, where m is the number of rows of A, that is, the  length of x. Therefore, each iteration on j will cost O(m) operations and this shows that the

the complexity of each iteration in k is O(n*m).

Remark 2: I did not implemented this algorithm and, therefore, I do not know if it will really work. Nevertheless, I think that its underlying Idea could be explored.

@Dario Fasino -: Well I actually tried your approach a long time ago. In fact I tried RQ factorization. The approach does not give me a correct answer because, note that matrix R is a rectangular matrix of dimensions (m * n), where m > n. All though it is upper triangular, there are entire columns in the matrix that are filled with zero.

Sp while solving R p = b via forward substitution, there is no unique answer because R  is not full column rank.

All though the method gives me an answer, it turns to be incorrect, but somewhat close. to the correct answer.

@André Camargo -: Thank you to for the reply. Will try your approach.

@Rittich - I tried CG, that also takes takes a long time.

Since there are only 2% elements are non-zeros, how are you storing matrix A? When you take Cholesky factorization , have you accounted for fill-inns? Otherwise computationally it is not difficult problem as D and M changes with iterations.

@RC Mittal-: I store Matrix A (0.2 % fill) as a sparse matrix (well I tried both Eigen C++ and MAtrix toolkit java). I guess they are storing it column wise. However, when I compute the cholesky factorization of ADA^T (2% fill), I allocate the full dense storage matrix of dimensions of ADA^T. It may be interesting to note that, on computing the Cholesky factorization of ADA^T, the factorization has nearly 98% non zeros. 

@Dario-: Sorry for the long break between my replies. Upon reading your solution, I have few questions. Kindly correct me if I am wrong. Now since Q is no longer a square matrix, either QQ^T or Q^TQ is no longer going to be an identity matrix. My understanding of your method is that you make use of this fact in your method. So will it not fail?

I think you have lost the sparsity when you take Cholesky  factorization. I believe without factorization you apply any conjugate gradient type of method , should solve it. If need be apply a preconditioner.

@Mittal -: I did try the Conjugate gradient solve with diagonal preconditioning. The link to the method is provided in http://eigen.tuxfamily.org/dox/classEigen_1_1ConjugateGradient.html

It's perfomance was no better than Cholesky.

but then you have computational time advantage.

@Mittal- : I will retry the method, I have some reservations with my testing of CG method. But during my earlier testing, when I mentioned performance, I actually meant computational time. 

Never form A*D*A^T or its Cholesky decomposition explicitly. I think an iterative method is the way to go.

Does your problem occur from a (weighted) least-square problem? In that case it might be better to apply a specially tailored iterative method directly to A (and incorporate D).

For D=identity, the Golub-Kahan bidiagonalization method (https://en.wikipedia.org/wiki/Bidiagonalization) is a state-of-the-art approach for this kind of problems. So since D is diagonal (and hopefully nonsingular), you may apply it to A*sqrt(D).