OPR-computation-related linear algebra problem

[Ether]

Just installed Octave3.6.4_gcc4.6.2_20130408 on this computer.

Results:

GNU Octave, version 3.6.4

Octave was configured for "i686-pc-mingw32".

octave:1> tic; d = load('k:/MyOctave/d.dat'); toc
Elapsed time is 0.0625 seconds

octave:2> tic; N = load('k:/MyOctave/N.dat'); toc
Elapsed time is 90.0469 seconds

octave:3> tic; r = N \ d; toc
Elapsed time is 1.03125 seconds

octave:4> tic; r = N \ d; toc
Elapsed time is 0.953125 seconds

octave:5> tic; Ns = sparse(N); toc
Elapsed time is 0.0625 seconds
octave:6> tic; rs = Ns \ d; toc
Elapsed time is 0.1875 seconds

octave:7> tic; Ns = sparse(N); toc
Elapsed time is 0.046875 seconds
octave:8> tic; rs = Ns \ d; toc
Elapsed time is 0.171875 seconds

Anybody know why Octave takes so long to load N ?

Load N times, all on the same computer:

Delphi....0.6 seconds

C.........1.3 seconds

SciLab....1.7 seconds

RLab......3.7 seconds

Python....6.5 seconds

Octave...90.0 seconds

[Greg McKaskle]

LV times were on this computer

Win 7 Professional 64-bit
Xeon CPU E5-1620 3.6GHz
16G RAM

Greg McKaskle

[Ether]

Is anybody else running Octave on a machine with 32-bit-XP Pro?

Are you having the same 30 second delay for Octave to load, and 90 seconds to load the 12MB N.dat file?

[Nikhil Bajaj]

I can check after work in the lab this evening if no one gets to it first, Ether!

[Ether]

Originally Posted by Nikhil Bajaj

I can check after work in the lab this evening if no one gets to it first, Ether!

Thank you.

[Foster]

Welcome to RLaB. New users type `help INTRO'
RLaB version 2.1.05 Copyright (C) 1992-97 Ian Searle
RLaB comes with ABSOLUTELY NO WARRANTY; for details type `help warranty'
This is free software, and you are welcome to redistribute it under
certain conditions; type `help conditions' for details
> tic();
> d=readm("d.dat");
> toc()
     21.9
> tic();
> N=readm("N.dat");
> toc()
     28.8
> tic();
> x = solve(N,d,"S");
> toc()
     34.8
>

And right, if you do toc(); you don't get the output. Which isn't wasn't what I expected.

Xeon 2.93 Ghz cluster. But this code may not be doing any multithreading.

[Greg McKaskle]

CD let me attach again. I attached the things I intended for the previous post.

As with many of the analysis functions, this calls into a DLL, to the function InvMatrixLUDri_head. So it seems to be using LU decomposition. I think the matrix qualifies as sparse, so that helps with performance.

The direct solver was almost ten seconds.

Greg McKaskle

[Michael Hill]

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

[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 ""

[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?

[MikeE]

Sorry I'm a bit late to the party.

I'm running Octave 3.2.4, so an older version than flameout
Hardware is Dell E6420 laptop (CPU i7-2640M @ 2.8GHz) Win 7 64bit

octave:2> N = load('N.dat');
octave:3> d = load('d.dat');
octave:4> tic; Ns = sparse(N); toc
Elapsed time is 0.136 seconds.
octave:5> tic; r = Ns \ d; toc
Elapsed time is 0.096 seconds.
octave:6> tic; r = N \ d; toc
Elapsed time is 0.921 seconds.

So solving the sparse matrix is around 230ms and direct solution is around 900ms.

[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.

[Ether]

RLaB is not a contender for fastest speed, but it's definitely the tiniest linear algebra app out there. It weighs in at under 1.5MB.

Makes a wonderful pop-up for that quick calculation, or for high-school students just learning linear algebra.

Very easy to use... and free.

Welcome to RLaB. New users type `help INTRO'
RLaB version 2.1.05 Copyright (C) 1992-97 Ian Searle

This is free software, and you are welcome to redistribute it under
certain conditions; type `help conditions' for details

> tic(); N=readm("N.dat"); toc()
3.66

> d=readm("d.dat");

> tic(); x=N\d; toc()
28.5

> tic(); x=solve(N,d,"S"); toc()
14.5
>

The second solution method (x=solve(N,d,"S")) tells RLaB that the matrix is symmetric, so it uses LAPACK subroutine DSYTRF (Bunch-Kaufman diagonal pivoting) to solve, which cuts the time in half.