使用模数的Numpy矩阵幂/指数？

Question

能否使用numpy的linalg.matrix_power函数并加上模数，以使得元素不会增长到超过某个特定值？

- John Smith

---

你能定义一下模数的意思吗？ - Benjamin

取模运算等同于求余数。比如10 mod 3 = 1, 24 mod 5 = 4等。linalg.matrix_power虽然快，但我希望在元素变得过大之前能够对它们应用模运算。 - John Smith

啊，模数运算：http://zh.wikipedia.org/wiki/模除 - Benjamin

正确，但要与矩阵指数一起使用，以避免元素爆炸。 - John Smith

“模数”通常指复数的绝对值（而“模除”确实用于整数除法的余数）。 - Eric O. Lebigot

是的，modulus 是一个名词（某物的模：向量、复数等，“符号乘以模数格式”），而 modulo（常缩写为“mod”）则是另一种不同的东西（副词、分词？）。这将帮助我永远不再把余数称为模数，但我仍然找不到 NumPy 中矢量化的、逐元素幂运算的函数，就像内置的 pow(x, y, z=None, /) 一样，它相当于 x**y（有两个参数）或 x**y % z（有三个参数）。 - Tomasz Gandor

Answer 1

为了防止溢出，你可以利用这样一个事实：如果你先对每个输入数字取模，那么你将得到相同的结果；事实上：

为了防止溢出，你可以利用这样一个事实：如果你先对每个输入数字取模，那么你将得到相同的结果；事实上：

(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

和

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

因此，您可以将幂次方 k 分解为 a+b+… 或 a*b*… 或两者的任意组合。上述等式使您仅需计算小数之间的幂运算，从而大大降低了整数溢出的风险。

- Eric O Lebigot

---

除了模反元素，例如pow(a, -1, p)的工作方式不起作用。 - Gregory Morse

Answer 2

使用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

- Unapiedra

---

感谢您提到Numpy的错误，这个错误似乎也出现在np.dot中，而且持续了近十年(1.24.2)，我刚刚花了几个小时在基于Numpy的“平方幂”的实现中寻找一个bug，转换为dtype=object后可以得到正确的结果。 - flonk

Answer 3

对我有效的是在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

- Rexcirus

Answer 4

"

明显的方法有什么问题？

例如：

"

import numpy as np

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

- Joe Kington

---

也许原帖中使用了大指数，导致数据溢出问题。例如，在加密算法中，指数运算与取模结合通常用于大整数运算。 - wim

不适用于模反元，不能有效地处理大的子结果。 - Gregory Morse

Answer 5

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

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

- Tamas Hegedus