Simple Matrix class

[Paul McGuire]

I've posted a simple Matrix class on my website as a small-footprint
package for doing basic calculations on matrices up to about 10x10 in
size (no theoretical limit, but performance on inverse is exponential).
Includes:
- trace
- transpose
- conjugate
- determinant
- inverse
- eigenvectors/values (for symmetric matrices)
- addition and multiplication (with constant or other matrix)

Matrices are easily built from formatted strings, as in:

m = Matrix( """1 2 3 4
5 11 20 3
2 7 11 1
0 5 3 1""")

Pretty much a no-strings-attached license, just don't hassle me about
little things like warranty, support, merchantability , accuracy, etc.
See it at
http://www.geocities.com/ptmcg/pytho...html#matrix_py

-- Paul

[Casey Hawthorne]

Do you calcalate the matrix inversion, when you don't need to?

"Paul McGuire" <pt***@austin.r r.comwrote:

>I've posted a simple Matrix class on my website as a small-footprint
package for doing basic calculations on matrices up to about 10x10 in
size (no theoretical limit, but performance on inverse is exponential).
Includes:
- trace
- transpose
- conjugate
- determinant
- inverse
- eigenvectors/values (for symmetric matrices)
- addition and multiplication (with constant or other matrix)

Matrices are easily built from formatted strings, as in:

m = Matrix( """1 2 3 4
5 11 20 3
2 7 11 1
0 5 3 1""")

Pretty much a no-strings-attached license, just don't hassle me about
little things like warranty, support, merchantability , accuracy, etc.
See it at
http://www.geocities.com/ptmcg/pytho...html#matrix_py

-- Paul

--
Regards,
Casey

[Paul McGuire]

On Jan 23, 4:05 pm, Casey Hawthorne <caseyhHAMMER_T [e-mail]>
wrote:

Do you calcalate the matrix inversion, when you don't need to?

No, the inversion is only calculated on the first call to inverse(),
and memoized so that subsequent calls return the cached value
immediately. Since the Matrix class is mutable, the cache is
invalidated if any of the matrix elements are updated. So after
changing a matrix element, the inversion would be recalculated at the
next call to inverse().

Hope that clears things up. You can also do your own experiments,
including adding verbose logging to the memoizing decorators - an
example is included in the source code. Also the footprint is quite
small, just one Python source file, about 600-700 lines of code, and
all Python, so 100% portable.

-- Paul

[Robert Kern]

Paul McGuire wrote:

I've posted a simple Matrix class on my website as a small-footprint
package for doing basic calculations on matrices up to about 10x10 in
size (no theoretical limit, but performance on inverse is exponential).

Why is that? A simple and robust LU decomposition should be no more than O(n**3).

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

[Paul McGuire]

On Jan 23, 6:59 pm, Robert Kern <robert.k...@gm ail.comwrote:

Paul McGuire wrote:

I've posted a simple Matrix class on my website as a small-footprint
package for doing basic calculations on matrices up to about 10x10 in
size (no theoretical limit, but performance on inverse is exponential).

Why is that? A simple and robust LU decomposition should be no more than O(n**3).

--
Robert Kern

Well "3" is an exponent isn't it? :)

In truth, in my laziness, I didn't *actually* test the performance.
But after measuring inversion times for nxn matrices for n=2 to 12, I
get these results (run on Python 2.4.2, on a 2GHz CPU):

n seconds ln(seconds)
2 0.0004112254490 45 -7.79636895604
3 0.0010224763203 1 -6.88552782893
4 0.0043754164286 2 -5.43175358002
5 0.0146999129778 -4.21991370509
6 0.0507813143849 -2.98022681913
7 0.143077961026 -1.94436561528
8 0.39962257773 -0.917234732978
9 1.14412558021 0.134640659841
10 3.01953516439 1.10510290046
11 8.76039971561 2.17024153354
12 21.8032182861 3.0820575867

Plotting n vs. ln(seconds) gives a fairly straight line of slope about
1.09, and exp(1.09) = 2.97, so your big-O estimate seems to line up
nicely with the experimental data - I couldn't have fudged it any
better.

-- Paul

[Gabriel Genellina]

At Tuesday 23/1/2007 22:33, Paul McGuire wrote:

>On Jan 23, 6:59 pm, Robert Kern <robert.k...@gm ail.comwrote:

Paul McGuire wrote:

I've posted a simple Matrix class on my website as a small-footprint
package for doing basic calculations on matrices up to about 10x10 in
size (no theoretical limit, but performance on inverse is exponential).

Why is that? A simple and robust LU decomposition should be no

more than O(n**3).

Well "3" is an exponent isn't it? :)

But constant!
x**2 is a "power" (or quadratic, or polynomial) function. 2**x is an
"exponentia l" function. They're quite different.

>In truth, in my laziness, I didn't *actually* test the performance.
But after measuring inversion times for nxn matrices for n=2 to 12, I
get these results (run on Python 2.4.2, on a 2GHz CPU):

n seconds ln(seconds)
2 0.0004112254490 45 -7.79636895604
3 0.0010224763203 1 -6.88552782893
4 0.0043754164286 2 -5.43175358002
5 0.0146999129778 -4.21991370509
6 0.0507813143849 -2.98022681913
7 0.143077961026 -1.94436561528
8 0.39962257773 -0.917234732978
9 1.14412558021 0.134640659841
10 3.01953516439 1.10510290046
11 8.76039971561 2.17024153354
12 21.8032182861 3.0820575867

Plotting n vs. ln(seconds) gives a fairly straight line of slope about
1.09, and exp(1.09) = 2.97, so your big-O estimate seems to line up
nicely with the experimental data - I couldn't have fudged it any
better.

Nope, such semilog plot shows that time grows exponentially, like
t=3*exp(n), and that's really bad! :(
The points should be aligned on a log-log plot to be a power function.
As Robert Kern stated before, this problem should be not worse than
O(n**3) - how have you implemented it?
--
Gabriel Genellina
Softlab SRL


_______________ _______________ _______________ _____
Preguntá. Respondé. Descubrí.
Todo lo que querías saber, y lo que ni imaginabas,
está en Yahoo! Respuestas (Beta).
¡Probalo ya!
http://www.yahoo.com.ar/respuestas

[Paul McGuire]

The points should be aligned on a log-log plot to be a power function.

As Robert Kern stated before, this problem should be not worse than
O(n**3) - how have you implemented it?

Sure enough, the complete equation is t = 5e-05exp(1.1n), or t = 5e-05
X 3**n.

As for the implementation, it's pretty much brute force. Here is the
code itself - the cacheValue decorator memoizes the calculated inverse
for the given Matrix.

@cacheValue
def det(self):
"Function to return the determinant of the matrix."
if self.isSquare() :
if self.numRows() 2:
multiplier = 1
firstRow = self[1]
tmp = self._rows[1:]
rangelentmp = range(len(tmp))
col = 0
detsum = 0
for val in firstRow:
if val:
#~ tmp2 = Matrix([
RowVector(t[0:col]+t[col+1:]) for t in tmp ])
tmp2 = self.getCachedM atrix([
RowVector(t[0:col]+t[col+1:]) for t in tmp ])
detsum += ( multiplier * val * tmp2.det() )
multiplier = -multiplier
col += 1
return detsum
if self.numRows() == 2:
return self[1][1]*self[2][2]-self[1][2]*self[2][1]
if self.numRows() == 1:
return self[1][1]
if self.numRows() == 0:
return 0
else:
raise MatrixException ("can only compute det for square
matrices")

def cofactor(self,i ,j):
i-=1
j-=1
#~ tmp = Matrix([ RowVector(r[:i]+r[i+1:]) for r in
(self._rows[:j]+self._rows[j+1:]) ])
tmp = self.getCachedM atrix([ RowVector(r[:i]+r[i+1:]) for r in
(self._rows[:j]+self._rows[j+1:]) ])
if (i+j)%2:
return -tmp.det()
else:
return tmp.det()
#~ return (-1) ** (i+j) * tmp.det()

@cacheValue
def inverse(self):
if self.isSquare() :
if self.det() != 0:
ret = Matrix( [ RowVector( [ self.cofactor(i ,j) for j
in self.colrange() ] )
for i in self.rowrange() ] )
ret *= (1.0/self.det())
return ret
else:
raise MatrixException ("cannot compute inverse for
singular matrices")
else:
raise MatrixException ("can only compute inverse for square
matrices")

[Gabriel Genellina]

At Wednesday 24/1/2007 02:40, Paul McGuire wrote:

The points should be aligned on a log-log plot to be a power function.
As Robert Kern stated before, this problem should be not worse than
O(n**3) - how have you implemented it?

Sure enough, the complete equation is t = 5e-05exp(1.1n), or t = 5e-05
X 3**n.

So either reimplement it better, or place a notice in BIG letters to
your users about the terrible time and memory penalties of using inverse()
--
Gabriel Genellina
Softlab SRL


_______________ _______________ _______________ _____
Preguntá. Respondé. Descubrí.
Todo lo que querías saber, y lo que ni imaginabas,
está en Yahoo! Respuestas (Beta).
¡Probalo ya!
http://www.yahoo.com.ar/respuestas

[Paul McGuire]

On Jan 24, 11:21 am, Gabriel Genellina <gagsl...@yahoo .com.arwrote:

At Wednesday 24/1/2007 02:40, Paul McGuire wrote:

The points should be aligned on a log-log plot to be a power function.
As Robert Kern stated before, this problem should be not worse than
O(n**3) - how have you implemented it?

Sure enough, the complete equation is t = 5e-05exp(1.1n), or t = 5e-05
X 3**n.

So either reimplement it better, or place a notice in BIG letters to
your users about the terrible time and memory penalties of using inverse()

Well, really, the "terrible time and memory penalties" aren't all that
terrible for matrices up to, oh, 10x10 or so, but you wouldn't want to
use this on anything much larger than that. Hey, guess what?! That's
exactly what I said in my original post!

And the purpose/motivation for "reimplemen ting it better" would be
what, exactly? So I can charge double for it? As it was, I felt a
little guilty for posting a solution to such a likely homework
assignment, but now there's room for a student to implement a kickass
inverter routine in place of the crappy-but-useful-for-small-matrices
brute force one I've written.

Cubs win!
-- Paul

[Robert Kern]

Paul McGuire wrote:

And the purpose/motivation for "reimplemen ting it better" would be
what, exactly? So I can charge double for it?

So you can have accurate results, and you get a good linear solver out of the
process. The method you use is bad in terms of accuracy as well as efficiency.

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco