4
$\begingroup$

I wrote the code below to invert an upper triangular matrix, avoiding any possible multiplication/subtraction by zero. It just uses $\frac{1}{6}n^3+\ldots$ flops instead of $n^3+\ldots$ flops.

function B = InvTrMat(A,n,type)
  % This code inverts an upper triangular matrix of order n
  % without doing any multiplication or subtraction by zero.
  B = zeros(n);

  for j = 1:n
    B(j,j) = 1 / A(j,j);
    for i = j-1:-1:1
      for k = i+1:j
        B(i,j) = B(i,j) - A(i,k) * B(k,j);
      end
      B(i,j) = B(i,j) * B(i,i);
    end
  end

  return
end

The problem is it is slower than the backslash command of Matlab. I used the following code to test and compare both:

n = 128;
t = zeros(1,100);
u = zeros(1,100);

for i = 1:100
  A = rand(n,n);
  A = A+A';
  A = A + n*eye(n);

  R = chol(A);

  tic;
  B = R \ eye(n);
  t(i) = toc;

  tic;
  B = InvTrMat(R,n,'u');
  u(i) = toc;
end

sum(t)/100
sum(u)/100

Does anyone know how to make my code faster than Matlab's code? It must be possible since I'm using just a fraction of flops from the other code.

My Matlab version is 7.10.0 (R2010a), installed on a Windows 8.0 64-bits PC.

Thank you, everyone.

$\endgroup$
  • $\begingroup$ Out of curiosity, what are the timings you get? I get about a factor of 83 difference. I'm actually surprised it's not more. $\endgroup$ – Doug Lipinski Nov 18 '14 at 1:33
11
$\begingroup$

MATLAB's \ (aka mldivide) command does not blindly compute the inverse of the matrix. Instead, it uses one of several algorithms based on the type of matrix (see the "Algorithms" section of http://www.mathworks.com/help/matlab/ref/mldivide.html). In the case of a triangular matrix, MATLAB will use a triangular solver which is at least as good as yours in terms of operation count (I haven't looked at your code too closely but they're probably the same algorithm).

The real difference is that the linear algebra routines in MATLAB are actually calls to highly optimized Fortran routines under the hood (i.e. LAPACK, see http://www.mathworks.com/company/newsletters/articles/matlab-incorporates-lapack.html). For any case that matches the intended use of these routines you will not be able to match this performance in native MATLAB code. This is due to the nature of the MATLAB language which receives only limited optimization on the fly by the just in time compiler.

For best performance in MATLAB, vectorize things and use the built in routines whenever possible.

$\endgroup$
  • $\begingroup$ The timings I get on a Intel i5-3210M laptop with 8GB of RAM are 2.6116e-004 for the backslash operator and 0.0135 for my code. But I'm glad to know the difference isn't that great in all computers. :) Thank you for your explanation. I'll bear in mind your comments when coding in MATLAB. $\endgroup$ – user11837 Nov 18 '14 at 15:08
  • $\begingroup$ Lapack itself is not highly optimized (though yes, being in fortran does help). It uses a lot of block algorithms that use level-3 BLAS, which helps, but another key ingredient is the BLAS, which is optimized. Matlab typically ships using the Intel MKL for its BLAS. $\endgroup$ – Stephen Mar 30 '15 at 3:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.