计算符号二进制矩阵的对数或逆矩阵的单个元素

Question

我正在尝试获取矩阵A 的一个修饰词A_adj 的单个元素，这两者都需要是符号表达式，其中符号x_i 是二进制的，矩阵A 是对称的和稀疏。 Python 的sympy 非常适合小问题：

from sympy import zeros, symbols

size = 4
A = zeros(size,size)
x_i = [x for x in symbols(f'x0:{size}')]

for i in range(size-1):
    A[i,i] += 0.5*x_i[i]
    A[i+1,i+1] += 0.5*x_i[i]
    A[i,i+1] = A[i+1,i] = -0.3*(i+1)*x_i[i]

A_adj_0 = A[1:,1:].det()

A_adj_0

这会计算辅因子矩阵的第一个元素A_adj_0（即对应的minor）并正确地给我0.125x_0x_1x_2 - 0.28x_2x_2^2 - 0.055x_1^2x_2 - 0.28x_1x_2^2，这是表达式我需要，但是有两个问题：

1. 这对于较大的矩阵是完全不可行的（对于 ~100 的 sizes，我需要这个）。
2. x_i 是二进制变量（即 0 或 1），sympy 似乎无法简化二进制变量的表达式，即简化多项式 x_i^n = x_i。

第一个问题可以通过求解线性方程组Ay = b 来部分解决，其中b 设置为第一个基向量[1, 0, 0, 0]，这样y 是@ 的倒数的第一列987654336@。 y 的第一个条目是A 的倒数的第一个元素：

b = zeros(size,1)
b[0] = 1

y = A.LUsolve(b)

s = {x_i[i]: 1 for i in range(size)}

print(y[0].subs(s) * A.subs(s).det())
print(A_adj_0.subs(s))

这里的问题是y的第一个元素的表达式非常复杂，即使在使用simplify()等之后也是如此。这将是一个非常简单的表达式，简化了上面第 2 点中提到的二进制表达式。这是一种更快的方法，但对于较大的矩阵仍然不可行。

这归结为我的实际问题：

是否有一种有效的方法来计算稀疏对称符号矩阵的对数的单个元素，其中符号是二进制值？

我也愿意使用其他软件。

附录 1：
似乎可以通过我不知道的简单自定义替换来简化sympy 中的二进制表达式：

A_subs = A_adj_0

for i in range(size):
    A_subs = A_subs.subs(x_i[i]*x_i[i], x_i[i])

A_subs

Answer 1

您应该确保在 sympy 中使用 Rational 而不是浮动，所以 S(1)/2 或 Rational(1, 2) 而不是 0.5。

sympy 中有一个新的（未记录的，目前是内部的）矩阵实现，称为 DomainMatrix。对于这样的问题，它可能要快得多，并且总是以完全扩展的形式产生多项式结果。我希望这种问题会更快，但它似乎仍然相当慢，因为它在内部并不稀疏（但是 - 这可能会在下一个版本中改变）并且它没有利用简化来自二进制值的符号。它可以在GF(2) 上工作，但不能使用假定在GF(2) 中的符号，这是不同的。

如果它有帮助，尽管这是您在 sympy 1.7.1 中使用它的方式：

from sympy import zeros, symbols, Rational
from sympy.polys.domainmatrix import DomainMatrix


size = 10
A = zeros(size,size)
x_i = [x for x in symbols(f'x0:{size}')]

for i in range(size-1):
    A[i,i] += Rational(1, 2)*x_i[i]
    A[i+1,i+1] += Rational(1, 2)*x_i[i]
    A[i,i+1] = A[i+1,i] = -Rational(3, 10)*(i+1)*x_i[i]

# Convert to DomainMatrix:
dM = DomainMatrix.from_list_sympy(size-1, size-1, A[1:, 1:].tolist())

# Compute determinant and convert back to normal sympy expression:
# Could also use dM.det().as_expr() although it might be slower
A_adj_0 = dM.charpoly()[-1].as_expr()

# Reduce powers:
A_adj_0 = A_adj_0.replace(lambda e: e.is_Pow, lambda e: e.args[0])

print(A_adj_0)