计算符号二进制矩阵的辅助矩阵或逆矩阵的单个元素
Computing a single element of the adjugate or inverse of a symbolic binary matrix
我正在尝试获取矩阵 A
的 adjugate 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,这是我需要的表达式,但是有两个问题:
- 这对于较大的矩阵是完全不可行的(我需要这个
size
s of ~100)。
x_i
是二进制变量(即0或1)并且sympy
似乎没有办法简化二进制变量的表达式,即简化多项式x_i^ n = x_i.
第一个问题可以通过求解线性方程组 Ay = b
来部分解决,其中 b
设置为第一个基向量 [1, 0, 0, 0]
,使得 y
是 A
的倒数的第一列。 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
你应该确保在 sympy 中使用 Rational 而不是浮动,所以 S(1)/2
或 Rational(1, 2)
而不是 0.5
。
在称为 DomainMatrix 的 sympy 中有一个新的(未记录且目前是内部的)矩阵实现。对于这样的问题,它可能要快得多,并且总是以完全展开的形式产生多项式结果。我希望它对于这类问题会快得多,但它似乎仍然相当慢,因为它在内部并不稀疏(但是 - 这可能会在下一个版本中改变)并且它没有利用简化来自二进制值的符号。它可以在 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)
我正在尝试获取矩阵 A
的 adjugate 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,这是我需要的表达式,但是有两个问题:
- 这对于较大的矩阵是完全不可行的(我需要这个
size
s of ~100)。 x_i
是二进制变量(即0或1)并且sympy
似乎没有办法简化二进制变量的表达式,即简化多项式x_i^ n = x_i.
第一个问题可以通过求解线性方程组 Ay = b
来部分解决,其中 b
设置为第一个基向量 [1, 0, 0, 0]
,使得 y
是 A
的倒数的第一列。 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
你应该确保在 sympy 中使用 Rational 而不是浮动,所以 S(1)/2
或 Rational(1, 2)
而不是 0.5
。
在称为 DomainMatrix 的 sympy 中有一个新的(未记录且目前是内部的)矩阵实现。对于这样的问题,它可能要快得多,并且总是以完全展开的形式产生多项式结果。我希望它对于这类问题会快得多,但它似乎仍然相当慢,因为它在内部并不稀疏(但是 - 这可能会在下一个版本中改变)并且它没有利用简化来自二进制值的符号。它可以在 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)