Path: csiph.com!fu-berlin.de!uni-berlin.de!not-for-mail
From: Peter Otten <__peter__@web.de>
Newsgroups: comp.lang.python
Subject: Re: LU decomposition
Date: Sun, 01 Nov 2015 12:55:50 +0100
Organization: None
Lines: 44
Message-ID: <mailman.19.1446378964.4463.python-list@python.org>
References: <1e6c4ce1-885c-43ac-b0a5-054f46d4c96e@googlegroups.com>
Mime-Version: 1.0
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 7Bit
User-Agent: KNode/4.13.3
Precedence: list
Xref: csiph.com comp.lang.python:98040

gers antifx wrote:

> I have to write a LU-decomposition. My Code worked so far but (I want to
> become better:) ) I want to ask you, if I could write this
> LU-decomposition in a better way?
> 
> def LU(x):
>     L = np.eye((x.shape[0]))
>     n = x.shape[0]
>     for ii in range(n-1):
>         for ll in range(1+ii,n):
>             factor = float(x[ll,ii])/x[ii,ii]
>             L[ll,ii] = factor
>             for kk in range(0+ii,n):
>                     x[ll,kk] = x[ll,kk] - faktor*x[ii,kk]
>     LU = np.dot(L,x)

You want to become better? A good way to find the problems Chris pointed out 
quickly is to write unit tests similar to

import unittest
...
from mymodule import LU

class LUTest(unittest.TestCase):
    def test_LU(self):
        matrix = ...
        expected_result = ...
        self.assertEqual(LU(matrix), expected_result)

if __name__ == "__main__":
   unittest.main()

There's no point making stylistic changes or tweaks to improve performance 
on code that doesn't work and doesn't have a well-defined interface yet,
but as a rule of thumb for efficient number-crunching you should try to 
replace for loops with numpy's array operations if at all possible.
Google found me 

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.lu.html

but it's hard to learn from that code as the real meat is a few levels down 
and implemented in C.