Numpy 矩阵幂/指数与模？

Question

是否可以将 numpy 的 linalg.matrix_power 与模一起使用，以便元素不会增长到大于某个值？

Answer 1

使用 Numpy 的实现：

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

我通过添加模项来对其进行调整。 但是，存在一个错误，即如果发生溢出，则不会引发 OverflowError 或任何其他类型的异常。从那时起，解决方案将是错误的。 此处有一个错误报告。

这是代码。小心使用：

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

Answer 2

这种显而易见的方法有什么问题吗？

例如

import numpy as np

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

Answer 3

我之前的所有解决方案都存在溢出问题，因此我必须编写一个算法来解决每次整数乘法后的溢出问题。我就是这样做的：

def matrix_power_mod(x, n, modulus):
    x = np.asanyarray(x)
    if len(x.shape) != 2:
        raise ValueError("input must be a matrix")
    if x.shape[0] != x.shape[1]:
        raise ValueError("input must be a square matrix")
    if not isinstance(n, int):
        raise ValueError("power must be an integer")

    if n < 0:
        x = np.linalg.inv(x)
        n = -n
    if n == 0:
        return np.identity(x.shape[0], dtype=x.dtype)
    y = None
    while n > 1:
        if n % 2 == 1:
            y = _matrix_mul_mod_opt(x, y, modulus=modulus)
        x = _matrix_mul_mod(x, x, modulus=modulus)
        n = n // 2
    return _matrix_mul_mod_opt(x, y, modulus=modulus)


def matrix_mul_mod(a, b, modulus):
    if len(a.shape) != 2:
        raise ValueError("input a must be a matrix")
    if len(b.shape) != 2:
        raise ValueError("input b must be a matrix")
    if a.shape[1] != a.shape[0]:
        raise ValueError("input a and b must have compatible shape for multiplication")
    return _matrix_mul_mod(a, b, modulus=modulus)


def _matrix_mul_mod_opt(a, b, modulus):
    if b is None:
        return a
    return _matrix_mul_mod(a, b, modulus=modulus)


def _matrix_mul_mod(a, b, modulus):
    r = np.zeros((a.shape[0], b.shape[1]), dtype=a.dtype)
    bT = b.T
    for rowindex in range(r.shape[0]):
        x = (a[rowindex, :] * bT) % modulus
        x = np.sum(x, 1) % modulus
        r[rowindex, :] = x
    return r

Answer 4

对我有用的是在一行中更改 Tamas Hegedus 代码，以指定 dtype=object。这是必需的，因为 Numpy 不像 python 那样支持大整数。使用 dtype=object 强制 numpy 支持大整数并避免无提示溢出错误。
使用指数 10^(15) 和模数约 10^(12) 进行测试。

def matrix_power_mod(x, n, modulus):
    x = np.asanyarray(x)
    if len(x.shape) != 2:
        raise ValueError("input must be a matrix")
    if x.shape[0] != x.shape[1]:
        raise ValueError("input must be a square matrix")
    if not isinstance(n, int):
        raise ValueError("power must be an integer")

    if n < 0:
        x = np.linalg.inv(x)
        n = -n
    if n == 0:
        return np.identity(x.shape[0], dtype=x.dtype)
    y = None
    while n > 1:
        if n % 2 == 1:
            y = _matrix_mul_mod_opt(x, y, modulus=modulus)
        x = _matrix_mul_mod(x, x, modulus=modulus)
        n = n // 2
    return _matrix_mul_mod_opt(x, y, modulus=modulus)


def matrix_mul_mod(a, b, modulus):
    if len(a.shape) != 2:
        raise ValueError("input a must be a matrix")
    if len(b.shape) != 2:
        raise ValueError("input b must be a matrix")
    if a.shape[1] != a.shape[0]:
        raise ValueError("input a and b must have compatible shape for multiplication")
    return _matrix_mul_mod(a, b, modulus=modulus)


def _matrix_mul_mod_opt(a, b, modulus):
    if b is None:
        return a
    return _matrix_mul_mod(a, b, modulus=modulus)


def _matrix_mul_mod(a, b, modulus):
    r = np.zeros((a.shape[0], b.shape[1]), dtype=object)
    bT = b.T
    for rowindex in range(r.shape[0]):
        x = (a[rowindex, :] * bT) % modulus
        x = np.sum(x, 1) % modulus
        r[rowindex, :] = x
    return r

Answer 5

为了防止溢出，您可以利用以下事实：如果您首先对每个输入数字取模，则会得到相同的结果；事实上：

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

对于一个矩阵 M。这来自以下两个基本恒等式，它们对于整数 x 和 y（以及正幂 p）有效：

(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*

相同的恒等式也适用于矩阵，因为矩阵加法和乘法可以通过标量加法和乘法来表达。这样，您只需对小数求幂（n mod p 通常比 n 小得多），并且发生溢出的可能性要小得多。因此，在 NumPy 中，您只需执行以下操作

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

即可获得 (arr**k) mod p。

如果这仍然不够（即，如果尽管 n mod p 很小，但仍存在 [n mod p]**k 导致溢出的风险），您可以中断将幂分解为多次幂。上述基本恒等式产生

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

and

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

因此，您可以将幂 k 分解为 a+b+… 或 a*b*… 或其任意组合。上述恒等式允许您仅执行小数乘小数的幂运算，这大大降低了整数溢出的风险。