Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #98031 > unrolled thread

LU decomposition

Back to article view | Back to comp.lang.python

Contents

  LU decomposition gers antifx <schweiger.gerald@gmail.com> - 2015-11-01 02:04 -0800
    Re: LU decomposition Chris Angelico <rosuav@gmail.com> - 2015-11-01 21:18 +1100
    Re: LU decomposition Peter Otten <__peter__@web.de> - 2015-11-01 12:55 +0100
    Re: LU decomposition Oscar Benjamin <oscar.j.benjamin@gmail.com> - 2015-11-03 15:52 +0000
    Re: LU decomposition Nagy László Zsolt <gandalf@shopzeus.com> - 2015-11-03 19:21 +0100

#98031 — LU decomposition

Hey,

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)

Thanks

[toc] | [next] | [standalone]

#98035

On Sun, Nov 1, 2015 at 9:04 PM, gers antifx <schweiger.gerald@gmail.com> 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)

For a start, I would recommend being careful with your variable names.
You have 'factor' all except one, where you have 'faktor' -
transcription error, or nearly-identical global and local names? And
all your other names are fairly opaque; can they be better named? I'm
particularly looking at this line:

x[ll,kk] = x[ll,kk] - faktor*x[ii,kk]

It is *extremely not obvious* that the first two are using 'll' and
the last one is using 'ii'. (Though I would write this as "x[ll,kk] -=
faktor*x[ii,kk]", which at least cuts down the duplication, so it's
less glitchy to have one out of three that's different.) I was going
to ask if you had some reason for not inverting the factor and simply
using "x[ll,kk] *= factor", till I read it again and saw the
difference.

I'm seeing here a lot of iteration over ranges, then subscripting with
that. Is it possible to iterate over the numpy array itself instead?
That's generally a more Pythonic way to do things.

Assigning to the local name LU at the end of the function seems odd.
Did you intend to return the dot-product?

Beyond that, I'd need a better comprehension of the mathematics behind
this, to evaluate what it's doing. So I'll let the actual experts dive
in :)

ChrisA

[toc] | [prev] | [next] | [standalone]

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

[toc] | [prev] | [next] | [standalone]

#98168

On 1 November 2015 at 10:04, gers antifx <schweiger.gerald@gmail.com> 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're using ii and ll, kk etc. This is common practise in Matlab code
because Matlab uses i as sqrt(-1) but this is spelled 1j in Python so
we don't have that problem. In this sort of context it would be more
natural to use simply i, j, k in Python.

Another is that you're not checking for divide by zero in
float(x[ll,ii])/x[ii,ii]. Are you sure that it will not be zero? What
happens if e.g. x[0, 0] is zero? This is where you'll need to use
pivoting which complicates the algorithm a bit.

--
Oscar

[toc] | [prev] | [next] | [standalone]

#98224

> Hey,
>
> 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?
Why can't you just use scipy for this?

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.lu.html

Do you really need to reinvent the wheel?

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web