OPR-computation-related linear algebra problem

[Nikhil Bajaj]

Originally Posted by Michael Hill

The reason I say it's computationally intensive is this article: http://www.johndcook.com/blog/2010/0...t-that-matrix/

That article is 100% correct. The solutions above that are solving in a handful of seconds or less are not inverting the matrix. Reducing the matrix to reduced-row echelon form is related to what methods like LU and Cholesky factorization do.

Even normal Gaussian elimination will be pretty fast on a sparse matrix (though still slower than most methods above), but it has problems with numerical stability that get worse and worse as matrix size increases and is for that main reason avoided by most people solving numerical linear algebra problems.

[DMetalKong]

Re-ran using Scipy's sparse matrix solver.

Average run time: 0.085s
Standard deviation: 0.005s

import sys
import numpy
import time
import scipy
import scipy.sparse
import scipy.sparse.linalg
import psutil

n_runs = 1000

print ""
print ""
print "Python version %s" % (sys.version)
print "Numpy version %s" % (numpy.__version__)
print "Scipy version %s" % (scipy.__version__)
print "Psutil version %s" % (psutil.__version__)
print ""


N = numpy.loadtxt(open('N.dat'))
d = numpy.loadtxt(open('d.dat'))

Ns = scipy.sparse.csr_matrix(N)

data = []
for i in range(1,n_runs+1):
    start = time.time()
    x = scipy.sparse.linalg.spsolve(Ns,d)
    end = time.time()
    row = [end - start]
    row.extend(psutil.cpu_percent(interval=1,percpu=True))
    s = "\t".join([str(item) for item in row])
    data.append(s)
    
f = open('times2.dat','w')
f.write("\n".join(data))
f.close()

_x = scipy.sparse.linalg.spsolve(Ns,d)
print ", ".join([str(f) for f in _x])
print ""

[James Critchley]

Nikhil,
Is there a reason why you are not using the "pcg" function which assumes symmetric positive definite inputs? This should be faster. Also please consider using the diagonal as a preconditioner. Unfortunately I do not have access the MATLAB at the moment. Could you try please the following? And sorry in advance for any bugs:

Ns = sparse(N);
D = diag(Ns);
Ds = sparse(diag(D)); #This was a bug... maybe it still is!

# Reference Solution
tic
output = Ns\d;
toc

# CG Solution
tic
output = pcg(Ns,d)
toc

# Diagonal PCG Solution
tic
output = pcg(Ns,d,[],[],Ds)
toc

# Reverse Cutthill-McKee re-ordering
tic
p = symrcm(Ns); # permutation array
Nr = Ns(p,p); # re-ordered problem
toc

# Re-ordered Solve
tic
output = Nr\d; #answer is stored in a permuted matrix indexed by 'p'
toc

Another advantage to the conjugate gradient methods is concurrent form of the solution within each iteration (parallel processing).

Best regards

[Nikhil Bajaj]

New Code, based on what James put up (I just added some disp's so that the results would be more clear. disps are outside of tics and tocs. I did not find any bugs though had to change #'s to %'s.

clc
disp('Loading Data...')
tic
d = load('d.dat');
N = load('N.dat');
toc
Ns = sparse(N);
D = diag(Ns);
Ds = sparse(diag(D)); %This was a bug... maybe it still is!

% Reference Solution 
disp('Reference Solution:')
tic
output1 = Ns\d;
toc


% CG Solution
disp('CG Solution:');
tic
output2 = pcg(Ns,d);
toc

% Diagonal PCG Solution
disp('Diagonal PCG Solution:');
tic
output3 = pcg(Ns,d,[],[],Ds);
toc

% Reverse Cutthill-McKee re-ordering
disp('Re-ordering (Reverse Cutthill-McKee:');
tic
p = symrcm(Ns); % permutation array
Nr = Ns(p,p); % re-ordered problem
toc

% Re-ordered Solve
disp('Re-ordered Solution:');
tic
output4 = Nr\d; %answer is stored in a permuted matrix indexed by 'p'
toc

Output:

Loading Data...
Elapsed time is 3.033846 seconds.
Reference Solution:
Elapsed time is 0.014136 seconds.
CG Solution:
pcg stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 4.8e-05.
Elapsed time is 0.007545 seconds.
Diagonal PCG Solution:
pcg converged at iteration 17 to a solution with relative residual 8.9e-07.
Elapsed time is 0.009216 seconds.
Re-ordering (Reverse Cutthill-McKee:
Elapsed time is 0.004523 seconds.
Re-ordered Solution:
Elapsed time is 0.015021 seconds.

I didn't precondition earlier because I was being sloppy/lazy . Thanks for calling me out. And you're right, I should have used pcg. Thanks for the suggestion.

[flameout]

Originally Posted by Ether

PS - can someone with a working Octave installation please run this? also SciLab and R

Since no-one has done Octave yet, I'll go ahead and do it (along with MATLAB for comparison). I can't do SciLab or R because I don't know how to use those

MATLAB 2012b:

>> N = dlmread('N.dat');
>> d = dlmread('d.dat');
>> tic ; r = N \ d; toc
Elapsed time is 0.797772 seconds.

GNU Octave 3.6.2:

octave:1> N = dlmread('N.dat');
octave:2> d = dlmread('d.dat');
octave:3> tic ; r = N \ d; toc
Elapsed time is 0.624047 seconds.

This is on an Intel i5 (2 core + hyperthreading) with Linux as the host OS (kernel version 3.7.6).

[Ether]

Originally Posted by flameout

Since no-one has done Octave yet, I'll go ahead and do it

Thanks flameout.

I can't do SciLab or R because I don't know how to use those

Just installed SciLab 5.4.1 with Intel Math Kernel Library 10.3 on a 7-year-old desktop PC:

• Intel Pentium D 3.4GHz (x86 Family 15 Model 6 Stepping 4)
  
• 32-bit XP Pro SP3
  
• 500GB Seagate Barracuda 7200

-->stacksize(70000000);
 
-->tic; N=read("N.dat",2509,2509); toc
 ans  = 1.672  
 
-->d=read("d.dat",2509,1);
 
-->tic; x=N\d; toc
 ans  = 1.672  
 
-->tic; Ns=sparse(N); toc
 ans  = 0.141  
 
-->tic(); xs = umfpack(Ns,'\',d); toc
 ans  = 0.14

[Ether]

Originally Posted by RyanCahoon

C code implementing Cholesky decomposition-based solver. With minimal optimization, the calculation runs in 3.02 seconds on my system.

Hi Ryan. I compiled it with Borland C++ 5.5 and ran it on the computer described in this post. It took 80 seconds:

ryan.exe 2509 N.dat d.dat x.dat

Reading: 1.312000 seconds
Calculation: 79.953000 seconds

So I dug up an old piece of code I wrote back in 1990 with a Cholesky factoring algorithm in it¹ and modified it for this application and ran it. It took about 22.5 seconds:

Nx=d build 5/26/2013 921p

CPU Hz (example 3.4e9 for 3.4GHz): 3.4e9
N matrix size (example 2509): 2509
N matrix filename (example N.dat): N.dat
d vector filename (example d.dat): d.dat
output filename (example x.dat): x.dat

reading N & d...
0.59 seconds

Cholesky...
22.37 seconds

Fwd & Back Subst...
0.08 seconds

Writing solution x...
0.01 seconds

done. press ENTER

If your code took only 3 seconds to run on your machine, but 80 on mine, I'm wondering what the Rice algorithm would do on your machine.


¹John Rischard Rice, Numerical Methods, Software, and Analysis, 1983, Page 139 (see attachments)

[flameout]

Originally Posted by Ether

Just installed SciLab 5.4.1 with Intel Math Kernel Library 10.3 on a 7-year-old desktop PC:

Now that I have working SciLab code, I'll go ahead and re-do the tests (with additional instrumentation on the reading). This is on the same computer as before (dual-core, hyperthreaded Intel i5, Linux 3.7.6).

MATLAB R2012b:

>> tic ; N = dlmread('N.dat'); toc
Elapsed time is 3.074810 seconds.
>> tic ; d = dlmread('d.dat'); toc
Elapsed time is 0.006744 seconds.
>> tic ; r = N \ d; toc
Elapsed time is 0.323021 seconds.
>> tic ; dlmwrite('out.dat', r); toc
Elapsed time is 0.124947 seconds.

I noticed that the solve time was very different from my previous run. This may be due to dynamic frequency scaling -- the results in this post (for all software) is with the frequency locked at the highest setting, 2.501 Ghz. It may also be due to a disk read -- I had not loaded MATLAB prior to running the previous test; now it's definitely in the disk cache. The solve is now consistently taking the time reported above, about a third of a second.

GNU Octave 3.6.2:

octave:1> tic ; N = dlmread('N.dat'); toc
Elapsed time is 1.87042 seconds.
octave:2> tic ; d = dlmread('d.dat'); toc
Elapsed time is 0.00241804 seconds.
octave:3> tic ; r = N \ d; toc
Elapsed time is 0.528489 seconds.
octave:4> tic ; dlmwrite('out.dat', r); toc
Elapsed time is 0.00820613 seconds.

Octave seems more consistent. The solve time is higher than for MATLAB, but the I/O times are consistently better.

Scilab 5.3.3:

-->stacksize(70000000);
 
-->tic; N=read("N.dat", 2509, 2509); toc
 ans  =
 
    1.21  
 
-->tic; d=read("d.dat", 2509, 1); toc
 ans  =
 
    0.003  
 
-->tic; x=N\d; toc
 ans  =
 
    1.052  
 
-->tic; Ns=sparse(N); toc
 ans  =
 
    0.081  
 
-->tic(); xs = umfpack(Ns,'\',d); toc
 ans  =
 
    0.081

Scilab failed to read the provided d.dat out of the box (reporting an EOF before it was done reading 2509 rows). I was able to correct this by adding a single newline to the end of d.dat.

FreeMat 4.0:

--> tic ; N = dlmread('N.dat'); toc
ans =
    2.8630
--> tic ; d = dlmread('d.dat'); toc
ans =
    0.0080
--> tic ; r = N \ d; toc
ans =
    3.4270

FreeMat did not have a dlmwrite function, so I haven't reported the write times for it. The time it took to solve the equations was significantly slower than any of the other programs. This did not improve with subsequent runs.

[Ether]

Originally Posted by Ether

Attached ZIP file contains N and d.

BTW, in case anyone was wondering...

Ax~b (overdetermined system)

A^TAx = A^Tb (normal equations; least squares solution of Ax~b)

Let N=A^TA and d=A^Tb (N is symmetric positive definite)

then Nx=d

A is the binary design matrix of alliances and b is the vector of alliance scores for all the qual matches for the 2013 season, including 75 Regionals and Districts, plus MAR and MSC, plus Archi, Curie, Galileo, and Newton.

So solving Nx=d for x is solving for 2013 World OPR.

[RyanCahoon]

Hi Ether,

Originally Posted by Ether

I compiled it with Borland C++ 5.5 and ran it [...] It took 80 seconds

That's quite a large difference in runtime. I compiled mine with Visual Studio 2010. I had wondered if VS was able to do any vectorized optimizations, but I don't see evidence of that in the Disassembly.

Originally Posted by Ether

If your code took only 3 seconds to run on your machine, but 80 on mine, I'm wondering what the Rice algorithm would do on your machine.

If I'm reading the pseudocode you posted correctly, I think I'm using the same algorithm (I got mine from the formulae on Wikipedia), the only difference I could find is I didn't handle the case of roundoff errors leading to slightly negative sums for the diagonal elements and I do some of the sums in reverse order, but unless there's some drastically bad cache effects I don't see that impacting the runtime.

Makes me wonder what you may have done better in your coding of the algorithm.

EDIT: Changing the order of the summations got me down to 2.68 and changing to in-place computation like your code got me to 2.58. Beyond that, any improvements would seem to be in the way the Pascal compiler is generating code.

Best,

[Ether]

Originally Posted by RyanCahoon

Makes me wonder what you may have done better in your coding of the algorithm.

...

[Greg McKaskle]

Finally had time to speak with the math guys.

The built-in LV linear algebra I was using links to an older version of Intel's MKL, but if I had used the SPD option on the solver it would indeed have been faster than the general version.

There is a toolkit called "Multicore Analysis and Sparse Matrix Toolkit", and they ran the numbers using that tool as well. Due to a newer version of MKL, the general solver is much faster. The right column converts the matrix into sparse form and uses a sparse solver.

Greg McKaskle

[Ether]

Originally Posted by Greg McKaskle

Finally had time to speak with the math guys...

Thanks Greg. Are the "time" units in the attachment milliseconds?

[Greg McKaskle]

Yes, they are in milliseconds. SPD stands for symmetric positive definite, column three enables the algorithms to utilize more than one core -- though this doesn't seem to help that much.

Greg McKaskle

[Ether]

I did the computation on this computer using a slightly modified version of DMetalKong's Python code.

Python 2.7.5
SciPy 0.12.0
NumPy 1.7.1

>>> import numpy
>>> import time
>>> import scipy
>>> import scipy.sparse
>>> import scipy.sparse.linalg
>>>
>>> # Read N & d ...
... start = time.time()
>>> N = numpy.loadtxt(open('E:\z\N.dat'))
>>> d = numpy.loadtxt(open('E:\z\d.dat'))
>>> end = time.time()
>>> print "%f seconds" % (end-start)
6.532000 seconds
>>>
>>> # solve...
... start = time.time()
>>> x = numpy.linalg.solve(N,d)
>>> end = time.time()
>>> print "%f seconds" % (end-start)
15.281000 seconds
>>>
>>> # Convert to sparse...
... start = time.time()
>>> Ns = scipy.sparse.csr_matrix(N)
>>> end = time.time()
>>> print "%f seconds" % (end-start)
0.234000 seconds
>>>
>>> # solve sparse...
... start = time.time()
>>> xs = scipy.sparse.linalg.spsolve(Ns,d)
>>> end = time.time()
>>> print "%f seconds" % (end-start)
0.453000 seconds
>>>

I had expected Python to be at least as fast as SciLab.

Perhaps there's an MKL for Python I need to install?