Numpy 2D Array - Power Of - Not returning an answer?

Question

I'm attempting to get the 'power' of a Python list/matrix using numpy. My only current working solution is an iterative function using np.dot():

def matr_power(matrix, power):
    matrix_a = list(matrix)
    matrix_b = list(matrix)
    for i in range(0, power-1):
        matrix_a = np.dot(matrix_a, matrix_b)
    return matrix_a

This works for my needs, but I'm aware it's probably not the most efficient method.

I've tried converting my list to a numpy array, performing power operations on it, and then back to a list so it's usable in the form I need. The conversions seem to happen, but the power calculation does not.

while (foo != bar):
    matr_x = np.asarray(matr_a)
    matr_y = matr_x ** n
    matr_out = matr_y.tolist()
    n += 1
    # Other code here to output certain results

The issue is, the matrix gets converted to an array as expected, but when performing the power operation (**) matr_y ends up being the same as matr_x as if no calculation was ever performed. I have tried using np.power(matr_y, n) and some other solutions found in related questions on Stack Overflow.

I've tried using the numpy documentation, but (either I'm misunderstanding it, or) it just confirms that this should be working as expected.

When checking the debugging console in PyCharm everything seems fine (all matrices / lists / arrays are converted as expected) except that the calculation matr_x ** i never seems to be calculated (or else never stored in matr_y).

Answer

Although it's possible to use a numpy matrix with the ** operator, the best solution is to use numpy arrays (as numpy matrices are deprecated) combined with numpy's linalg matrix_power method.

matr_x = np.array(mat_a)
matr_y = np.linalg.matrix_power(matr_x, path_length)
work_matr = matr_y.tolist()

It is also now apparent that the function of ** being element-wise may have been discovered earlier had I not been using an adjacency matrix (only zeros and ones).

Okay thanks, but power(matr_y, i) also doesn't seem to act as I expected, is this also element-wise? Is iterating through np.dot() my only option? — Matthew, Feb 01 '17 at 13:11
Yes, `power` is also element-wise, as explained in the first line of its docstring. :) — Warren Weckesser, Feb 01 '17 at 13:11

score 1 · Accepted Answer · answered Feb 01 '17 at 13:16

1

There are (at least) two options for computing the power of a matrix using numpy without multiple calls to dot:

Use numpy.linalg.matrix_power.
Use the numpy matrix class, which defines ** to be the matrix algebraic power.

For example,

In [38]: a
Out[38]: 
array([[0, 1, 0],
       [1, 0, 1],
       [0, 1, 0]])

In [39]: np.linalg.matrix_power(a, 2)
Out[39]: 
array([[1, 0, 1],
       [0, 2, 0],
       [1, 0, 1]])

In [40]: np.linalg.matrix_power(a, 3)
Out[40]: 
array([[0, 2, 0],
       [2, 0, 2],
       [0, 2, 0]])

In [41]: m = np.matrix(a)

In [42]: m ** 2
Out[42]: 
matrix([[1, 0, 1],
        [0, 2, 0],
        [1, 0, 1]])

In [43]: m ** 3
Out[43]: 
matrix([[0, 2, 0],
        [2, 0, 2],
        [0, 2, 0]])

answered Feb 01 '17 at 13:16

Warren Weckesser

110,654
19
194
214

1

Matthew, if you are unsure which one to use. Here is kind of the [official position](http://stackoverflow.com/a/20964252/7207392) on this. – Paul Panzer Feb 01 '17 at 13:27
1

@Matthew I was just linking that in case you were torn between the two perfectly fine solutions Warren has offered. Personally I tend to agree that the `matrix` class doesn't offer enough to offset the confusion it causes. So I would go with `linalg.matrix_power`. If you ever have to program an efficient integer power routine yourself: You can save a lot of multiplications by repeatedly squaring your operand (M, M^2, M^4, M^8 ...) and then multiplying those whose bit in the binary representation of the exponent is set. – Paul Panzer Feb 01 '17 at 13:50
@PaulPanzer that sounds really interesting, but I'm a little confused as to how that works. Any chance you can expand or link to a concise explanation? – Matthew Feb 01 '17 at 13:57

score 1 · Answer 2 · answered Feb 01 '17 at 14:48

Warren's answer is perfectly good.

Upon special request by the OP I briefly explain how to build an efficient integer power operator by hand.

I don't know what this algorithm is called, but it works like this: Suppose you want to calculate X^35. If you do that naively it will cost you 34 multiplications. But you can do much better than that. Write X^35 = X^32 x X^2 x X. What you've done here is split the product according to the binary representation of 35, which is 100011. Now, calculating X^32 is actually cheap, because you only have to repeatedly (5 times) square X to get there. So in total you need just 7 multiplications, much better than 34.

In code:

def my_power(x, n):
    out = None
    p = x
    while True:
        if n % 2 == 1:
            if out is None:
                out = p
            else:
                out = out @ p # this requires a fairly up-to-date python
                              # if yours is too old use np.dot instead
            if n == 1:
                return out
        n //= 2
        p = p @ p

Numpy 2D Array - Power Of - Not returning an answer?

Answer

2 Answers2