Fast Quadruple Precision Computations in MATLAB

by Pavel Holoborodko on January 20, 2013

In order to improve performance we have made important changes to the toolbox recently.

Now toolbox is heavily optimized for computations with quadruple precision.

In more details, new core engine of toolbox (written in Assembler/C++) has two different variants for every routine:

  • The first variant is truly precision independent, capable of computing results with any required precision (same as before).
  • The other is highly optimized for quadruple precision (128-bit or 34 decimal digits) computations.

Toolbox automatically switches to the fast quadruple routines if the user asks toolbox to operate with the precision of 34 decimal digits. No additional actions are required.

Quadruple precision routines have an advantage of knowing size of numbers in advance thus avoiding expensive memory allocations. Also this allows more efficient usage of hardware capabilities of modern CPUs, even exploiting some of the advanced instructions, not to mention possibilities of harnessing power of the GPU!

Quadruple precision numbers and operations are fully compliant with IEEE 754-2008 standard.

Comparison with Symbolic Math Toolbox

We compare optimized quadruple precision mode in our toolbox with Symbolic Math Toolbox from MathWorks. We use 64-bit MATLAB 2012b, Windows 7 and PC based on Core 2 Quad CPU.

At first we do initialization for both toolboxes to operate with quadruple precision:

% Symbolic Math Toolbox uses 9 guard digits. 
% Hence we should use 25 digits of precision 
% to get 'true' quadruple precision (34 = 25+9)
>> digits(25);           
 
% Setup 34-digits mode in Multiprecision Computing Toolbox.
>> mp.Digits(34);
>> mp.GuardDigits(0); 
 
% Show all digits of the numbers
>> format longG
 
% Prepare test matrices
>> A = rand(100);       % 100 x 100 full rank test matrix
>> X = mp(A);           
>> Y = vpa(A);

Then we do tests for svd, eig and qr functions (timing is highlighted with yellow):

Singular Value Decomposition

>> tic; [U,S,V] = svd(Y); toc;      % Symbolic Math Toolbox
Elapsed time is 59.357723 seconds. 
>> norm(Y-U*S*V',1)
3.061145275981987186191661e-31
 
>> tic; [U,S,V] = svd(X); toc;       % Multiprecision Computing Toolbox
Elapsed time is 0.668525 seconds.    % x80 times faster 
>> norm(X-U*S*V',1)                 
5.718086004885703247723198511534649e-31

Eigenvalues and Eigenvectors

>> tic; [V,D] = eig(Y); toc;        % Symbolic Math Toolbox
Elapsed time is 108.474148 seconds. 
>> norm(Y*V - V*D,1)
1.664003471950571777040365e-30
 
>> tic; [V,D] = eig(X); toc;         % Multiprecision Computing Toolbox
Elapsed time is 1.433337 seconds.    % x75 times faster 
>> norm(X*V - V*D,1)
3.520599938339867139386327679045247e-31

QR Decomposition

>> tic; [Q,R] = vpa_qr(Y); toc;      % Symbolic Math Toolbox
Elapsed time is 7.269699 seconds. 
>> norm(Y-Q*R,1)
3.608711233295583179663639e-31
 
>> tic; [Q,R] = qr(X); toc;          % Multiprecision Computing Toolbox
Elapsed time is 0.212965 seconds.    % x34 times faster 
>> norm(X-Q*R,1)
4.758852521959335590660505497505944e-32

* Function vpa_qr is borrowed from here.

Large Matrices Computing

Fast ‘quadruple’ precision mode even makes it possible to handle sufficiently large matrices in reasonable time:

>> A = mp(rand(10000,400));               % 10000 x 400 test matrix
 
>> tic; [U,S,V] = svd(A,'econ'); toc;
Elapsed time is 332.870674 seconds.       % 5.5 minutes
 
>> norm(A - U*S*V',1)                     % Check accuracy of the result
8.332401597108828366965178092631891e-29

By our estimation, Symbolic Math Toolbox requires 7-8 hours to compute the same result (in fact, we stopped it after 5 hours of computations).

Leave a Comment

Previous post:

Next post: