Matrix-matrix and matrix-vector multiplication with triangular matrices

[afossa]

The matrix-matrix product and matrix-vector product with triangular matrix obtained via Cholesky decomposition does not return the expected result due to the upper triangular part not being set equal to zero.

Considering a positive definite square matrix A and its Cholesky decomposition A = S*S', I expect the following snippet to return A2 equal to A:

Matrix S = A.uplo(UPLO.LOWER).cholesky().lu;
Matrix A2 = S.mt(S);

But in fact what I get back is a different matrix, presumably due to the fact that the upper triangular part of S is not zero. The same unexpected behaviour happens with matrix-vector multiplication:

double[] y = S.mv(x)

does not return the expected value for y.

Up to now, I came out with the solution of manually setting equal to zero all coefficients in the upper triangular part of S, so that if A is 3x3 the following works:

S.set(0, 1, 0.0).set(0, 2, 0.0).set(1, 2, 0.0).uplo(null);
Matrix A3 = S.mt(S)

and A3 is equal to A. Of course, this is not convenient, and I wonder if I am misunderstanding the usage of UPLO or if it is a bug in the implementation.

Thank you!

[haifengl]

It is LAPACK convention that the only UPLO parts contains the factorization and the other half is untouched.

[afossa]

Thank you for the answer. And there exist any way (excluding loops over matrix elements) to achieve one of the followings in smile?

1. set equal to zero the untouched part of the matrix
2. extract the lower triangular elements of a matrix into another one
3. perform matrix-vector and matrix-matrix multiplication with matrices resulting from Cholesky decomposition (with half of the matrix untouched) and obtain the result you would expect if the untouched part is set to zero

I tried different combinations of UPLO and Matrix.mm(), Matrix.mv() and Matrix.mt() so far, but without success. Thank you for your advice,

[haifengl]

Check out http://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html

[afossa]

I checked out BLAS documentation, and what I am looking for is trmv() to perform matrix-vector multiplication when the matrix is considered triangular and only one UPLO part is considered during computation.

However, for Matrix objects, even if declared triangular (uplo and diag not null), trmv() is used only in one specific case, namely Matrix.mv(Transpose trans, double alpha, double[] x, double beta, double[] y) with alpha equal to 1, beta equal to 0 and x == y (lines 1479-1493 in Matrix.java).

Is there any specific reason behind this choice? I would expect trmv() being called with any alpha, beta and x != y, so as to correctly perform multiplications with triangular matrices.

I thought about calling BLAS.engine.trmv(...) directly in my code, but this is not convenient because there is no direct public access to the matrix storage A.

Thanks!

[haifengl]

It is the only use case that trmv supports.