1

I have a Java program making linear algebra. I use the jblas library that, according to my understanding, is supposed to call native libraries implementing Blas and Lapack for faster results.

This code is running within a Docker and launched in AWS Batch managed jobs.

Excerpt from Dockerfile :

FROM debian:stretch
RUN apt-get -y update && apt-get -y upgrade
RUN apt-get -y install libgfortran3 # install fortran
RUN apt-get -y install openjdk-8-jdk # install java 8

I try to improve the speed to inverse a 24000 by 24000 square symmetric matrix. I see there is 2 methods provided by jblas library. One for general purpose linear system solving using native dgesv procedure and another one for symmetric matrix solving using native dsysv procedure.

DoubleMatrix A = ...; // my 24k symmetric square matrix
DoubleMatrix B = DoubleMatrix.eye(A.rows);// Identity matrix
DoubleMatrix C = Solve.solve(A, B);// takes 4020 s
DoubleMatrix D = Solve.solveSymmetric(A, B);// takes 5040 s, longer than the calculation of C

Question 1 : is it normal that solving a 24k square matrix takes so much time when native Blas and Lapack libraries are supposed to be used ? If no, how to improve the speed while running in the context of a AWS Batch job ?

Question 2 : Why solve symmetric (dsysv) is slower than general solve (dgesv) ? My expectation is that if we let the native library know that the input matrix is symmetric, it gives it a hint that should allow it to solve the linear system faster.

By the way, I checked that the 2 ways of doing gives the same numerical results. This is the case.

Comencau
  • 1,084
  • 15
  • 35

1 Answers1

1

I will answer only one of your questions.

  • is it normal that solving a 24k square matrix takes so much time. A: Yes, matrix inversion is O(n^3) so you are doing dozens of trillions of operations. What are you using the matrix inverse for? If you just want to solve a system then it may be better to use an iterative algorithms that resembles gradient descent methods popular for neural networks, e.g. Gauss-Seidel, or Richardson.

The other question is tricky, not all the operations in a symmetric matrix are more efficient than in a general matrix. For sake of the argument suppose that you want one single element in a general matrix in column-major 1D array form. A[i,j] = A_ge[j*N+i], but in the symmetric representation where the lower triangular part A[i,j] = i > j ? A_sy[i*(i+1)/2 + j] : A_sy[j*(j+1)/2 + i], the redundancy in the memory made it simpler to retrieve the desired element.

I think the answer should be a no if the most efficient algorithm was implemented. But the main motivation for having a symmetric matrix representation may be to save memory or numerical stability, and not speed.

This 20% increase in run-time may be due to less efficient indexing, for instance,

Bob
  • 13,867
  • 1
  • 5
  • 27