Potencia de matriz numpy/exponente con módulo?

Question

¿Es posible usar linalg.matrix_power de numpy con un módulo para que los elementos no crezcan más que un cierto valor?

Answer 1

Con el fin de evitar el desbordamiento, puede utilizar el hecho de que se obtiene el mismo resultado si se toma el primer módulo de cada uno de sus números de entrada; de hecho:

(M**k) mod p = ([M mod p]**k) mod p, 

para una matriz M. Esto viene de las siguientes dos identidades fundamentales, que son válidos para los números enteros x y y:

(x+y) mod p = ([x mod p]+[y mod p]) mod p # All additions can be done on numbers *modulo p* 
(x*y) mod p = ([x mod p]*[y mod p]) mod p # All multiplications can be done on numbers *modulo p* 

Las mismas identidades se mantienen para matrices así, puesto que la adición de matriz y la multiplicación se pueden expresar a través de la suma escalar y multiplicación. Con esto, solo expones números pequeños (n mod p generalmente es mucho más pequeño que n) y es mucho menos probable que tengan desbordamientos. En NumPy, tendría por lo tanto, simplemente no hacer

((arr % p)**k) % p 

el fin de obtener (arr**k) mod p.

Si esto aún no es suficiente (es decir, si existe el riesgo de que [n mod p]**k provoque un desbordamiento a pesar de que n mod p es pequeño), puede dividir la exponenciación en múltiples exponenciaciones. Las identidades fundamentales anteriores rendimiento

(n**[a+b]) mod p = ([{n mod p}**a mod p] * [{n mod p}**b mod p]) mod p 

y

(n**[a*b]) mod p = ([n mod p]**a mod p)**b mod p. 

Por lo tanto, se puede romper el poder k como a+b+… o a*b*… o cualquier combinación de los mismos. Las identidades anteriores le permiten realizar solo exponenciaciones de números pequeños por números pequeños, lo que reduce en gran medida el riesgo de desbordamientos de enteros.

Answer 2

¿Qué pasa con el enfoque obvio?

E.g.

import numpy as np 

x = np.arange(100).reshape(10,10) 
y = np.linalg.matrix_power(x, 2) % 50 

Answer 3

Uso de la aplicación de Numpy:

https://github.com/numpy/numpy/blob/master/numpy/matrixlib/defmatrix.py#L98

I adaptado mediante la adición de un término de módulo. SIN EMBARGO, hay un error, en el caso de que ocurra un desbordamiento, no se genera OverflowError ni ningún otro tipo de excepción. Desde ese momento, la solución será incorrecta. Hay un informe de error here.

Aquí está el código. Usar con cuidado:

from numpy.core.numeric import concatenate, isscalar, binary_repr, identity, asanyarray, dot 
from numpy.core.numerictypes import issubdtype  
def matrix_power(M, n, mod_val): 
    # Implementation shadows numpy's matrix_power, but with modulo included 
    M = asanyarray(M) 
    if len(M.shape) != 2 or M.shape[0] != M.shape[1]: 
     raise ValueError("input must be a square array") 
    if not issubdtype(type(n), int): 
     raise TypeError("exponent must be an integer") 

    from numpy.linalg import inv 

    if n==0: 
     M = M.copy() 
     M[:] = identity(M.shape[0]) 
     return M 
    elif n<0: 
     M = inv(M) 
     n *= -1 

    result = M % mod_val 
    if n <= 3: 
     for _ in range(n-1): 
      result = dot(result, M) % mod_val 
     return result 

    # binary decompositon to reduce the number of matrix 
    # multiplications for n > 3 
    beta = binary_repr(n) 
    Z, q, t = M, 0, len(beta) 
    while beta[t-q-1] == '0': 
     Z = dot(Z, Z) % mod_val 
     q += 1 
    result = Z 
    for k in range(q+1, t): 
     Z = dot(Z, Z) % mod_val 
     if beta[t-k-1] == '1': 
      result = dot(result, Z) % mod_val 
    return result % mod_val