1

I have found a code that pivotize a square matrix for LU decomposition, but I can't understand some of them.

def pivotize(m):
    """Creates the pivoting matrix for m."""
    n = len(m)
    ID = [[float(i == j) for i in xrange(n)] for j in xrange(n)]
    for j in xrange(n):
        row = max(xrange(j, n), key=lambda i: abs(m[i][j]))
        if j != row:
            ID[j], ID[row] = ID[row], ID[j]
    return ID

First, isn't the line for ID simply the identity matrix? any advantage doing this?

Second, I can't really understand the line for row. I know lambda is used to define a function in the text, and it simply returns the value of M_ij once the value of i is supplied (and the value of j depends on the for loop), but what is i?

And isn't xrange similar to range? But what does it return here?

And when combined with the function max, what happens? I just don't know what are the things inside the max function that are being compared.

Sorry if this question sounds stupid. I'm pretty new to programming

Physicist
  • 2,848
  • 8
  • 33
  • 62

2 Answers2

3

First, isn't the line for ID simply the identity matrix?

Yes.

Second, I can't really understand the line for row....

See this for a discussion about the max/key/lambda interaction. To answer "what is i?", its the argument to the lambda function, i could equivalently be x for foo. (For clarity, yielding abs(m[x][j]) and abs(m[foo][j]), respectively).

And isn't xrange similar to range?

Yes. In Python 2 xrange returns a sequence object that lazily computes the next value, only when needed. See this for more info. When looped through in entirety, range and xrange will produce the same results, with different implementations.

But what does it return here?

Line 5's xrange(n) will return the integer values from 0 to (n-1) while Line 6's xrange(j, n) will return the integer values from j to (n-1).

EDIT

More on lambda:

Consider how you might take a given sequence of numbers, and double each of them. First define a function that doubles a number x and returns that value. Then map that function to each element of a sequence.

# Traditional
def double(x): return x*2
print map(double, [1,2,3,4])         # [2, 4, 6, 8]

You could have alternatively used an anonymous (lambda) function to do the same thing:

# Lambda
print map(lambda i: i*2, [1,2,3,4])  # [2, 4, 6, 8]

Note that everything is the same, except the definition of the double function is gone and the reference to the function in the call to map is replaced by the "inline" lambda function.

EDIT 2

And when combined with the function max, what happens?

This line row = max(xrange(j, n), key=lambda i: abs(m[i][j])) can be broken down as follows:

  • xrange(j,n) generates a sequence of integers from j (inclusive) to n (exclusive).
  • Each of those integers are then "passed" as arguments to the function in the key argument. In other words, they're each used as the i in the lambda function. The lambda function "returns" the absolute value of the ith row and jth column.[1]
  • The max function then finds the greatest value of these "lambda outputs" and sets row equal to it.

This could have alternatively been written as a max of a list comprehension:

row = max( [abs(m[i][j]) for i in xrange(j,n)] )

Or as Dan D. points out in his comment, written as a generator expression (without creating an intermediary list) as simply:

row = max( abs(m[i][j]) for i in xrange(j,n) )

Notes:

[1] Some assumptions here, but row-column is a standard way of expressing matrices.

Community
  • 1
  • 1
jedwards
  • 29,432
  • 3
  • 65
  • 92
  • I know 'i' is the argument to the lambda function, but 'where' is i? Normally after I define a function (say y=np.sin(t)), I will declare t=np.linspace(0,2*np.pi) as the arguments, but now I'm just confused as there aren't any other codes related to i... – Physicist Mar 05 '15 at 00:09
  • re: your comment edit -- `i` is no more made up than `x` is in the definition of the `double` function: `double(x):` – jedwards Mar 05 '15 at 00:17
  • 1
    @Physicist: `i` is just an arbitray name given to the argument that gets passed to the `key` function whose return value determines the ordering used by the `max()` function. In this context it will be an integer with a value from `j` to `n-1` inclusive because of the `xrange(j, n)`. – martineau Mar 05 '15 at 00:25
  • 1
    Note that as `max` can take any iterable, even a generator expression, you can write `row = max( abs(m[i][j]) for i in xrange(j,n) )` and save creating an actual `list` due to the list comprehension. – Dan D. Mar 05 '15 at 00:30
2

xrange is an Python2 construct that's used for handling range memory efficient. Beforerange actually created a list, and then for loop would run through it. xrange however is a generator, that means it spits out 1 value at a time when it's asked to, without creating a full list.

ID is actually an identity matrix. You're right there. That's a neat trick that boolean can be converted to float of value 1.0.

Then the snippet runs through all the remaining rows and finds the maximum value of that row in the original matrix row = max(xrange(j, n), key=lambda i: abs(m[i][j])). Notice that there's a 2nd nice trick in there, max can take any iterable object to operate on, including generators. lambda keyword in that row indicates what's known as "anonymous function".

In more details: Anonymous functions are functions that are not bound to an identifier. Lambda function in your snippet takes in 1 value, i and returns an absolute value at matrix position m[i][j]. The values that get sent as the input for the function are provided by the generator xrange(j, n).

max then takes the return values of the lambda function as the thing it actually compares. In example, in python3 it's not possible to compare 2 different types. I.e. comparing string > int produces: TypeError: unorderable types: str() > int(). However if we're sure that the list contains numbers, just in different format, you can do something like:

>>> l = ["1", "2", 3, 4]
>>> max(l, key=lambda x: int(x))
4
>>> min(l, key=lambda x: int(x))
'1' #type 'str'

This just shows, that the actual compared values are the returns of your key function, but actual produced values are your original input values.

Once it found the rows maximum value, it "rotates", by substitution ID[j], ID[row] = ID[row], ID[j], all the other rows, in the identity matrix, around it so that only the maximum value remains on the diagonal.

This helps prevent division by a too small a number in the next step of LU decomposition.

What you get in return, isn't the original matrix, rotated, but a matrix of 1.0s and 0.0s that is your transformation matrix that multiplied by the original one will result in the pivoted matrix.

This seems a very nicely written function that saves memory and helps performance in python. Hopefully I did this right.

ljetibo
  • 3,048
  • 19
  • 25