Numpy matrix power/esponente con modulo?

Question

È possibile utilizzare linalg.matrix_power di numpy con un modulo in modo che gli elementi non aumentino di un certo valore?

Answer 1

Al fine di evitare overflow, è possibile utilizzare il fatto che si ottiene lo stesso risultato, se si prende il primo modulo di ciascuno dei tuoi numeri di ingresso; infatti:

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

per un matriceM. Ciò deriva dalle seguenti due identità fondamentali, che sono validi per gli interi x e 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* 

Le stesse identità valgono per matrici pure, poiché oltre matrice e moltiplicazione possono essere espresse attraverso scalari e della moltiplicazione. Con questo, si esponentano solo numeri piccoli (n mod p è generalmente molto più piccolo di n) e sono molto meno probabilità di ottenere overflow. In NumPy, si dovrebbe pertanto fare semplicemente

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

al fine di ottenere (arr**k) mod p.

Se questo non è ancora sufficiente (vale a dire, se v'è il rischio che [n mod p]**k causa di overflow nonostante n mod p essendo piccolo), è possibile suddividere l'elevazione a potenza in più exponentiations. Le identità fondamentali di cui sopra resa

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

e

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

Così, si può rompere il potere k come a+b+… o a*b*… o qualsiasi loro combinazione. Le identità sopra consentono di eseguire solo l'esponenziazione di piccoli numeri con numeri piccoli, il che riduce notevolmente il rischio di overflow di interi.

Answer 2

Cosa c'è di sbagliato nell'approccio ovvio?

E.g.

import numpy as np 

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

Answer 3

Usando l'implementazione Numpy:

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

ho adattato con l'aggiunta di un termine modulo. TUTTO IL, c'è un bug, nel senso che se si verifica un overflow, non viene generato alcun OverflowError o nessun altro tipo di eccezione. Da quel momento in poi, la soluzione sarà sbagliata. C'è un bug report here.

Ecco il codice. Usare con cautela:

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