将 sympy 表达式转换为向量以找到线性独立的子集

Turning a sympy expression into a vector to find linearly independent subset

我有一个表达式列表,例如 4.0*x[0] + 5.0*x[10] + 1 = 0 我想根据 [4.0, 0, 0, ..., 5.0, ... , 1] 等系数将它们转换为向量。原因是我的一些方程可能是线性相关的,我想从 numpy 库中 运行 QR 这样我就可以找到一个线性独立的子集。

我可以通过 expr.replace(x[i], 0)i 通配符索引来获得常数项。我还可以通过 expr.atoms(Mul) 得到大部分其他术语,这给了我集合 4.0*x[0], 5.0*x[10] 然后对于这些表达式中的每一个,我可以做 expr.atoms(Indexed).pop()expr.atoms(Float).pop() 来拆分部分.

问题是当我有一个像 x[0] + 5.0*x[10] + 1 = 0 这样的表达式时,其中第一个变量出现的隐式系数为 1。该术语不再被识别为 Mul 对象。

无论如何,我认为可能有更好的方法来实现我的目标?

如果你给你的符号一个特定的顺序,就像下面的代码一样,你可以将表达式转换为多项式并得到它的系数:

>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> a1, a2, a3, a4 = symbols('a[1], a[2], a[3], a[4]')
>>> used_symbols = (a1, a2, a3, a4)
>>> replacements = [(n, x**(enu+1)) for enu,n in enumerate(used_symbols)]
>>> expr = 5 + a1 + 4*a4
>>> Poly(expr.subs(replacements)).all_coeffs()
[4, 0, 0, 1, 5]

如果事先不知道使用以下递归函数,您也可以检索所用符号的列表:

def retrieve_used_symbols(expr):
    """Return the symbols used in the `expr` in a list."""
    used_symbols = []
    for term in expr.args:
        if term.is_Atom and term.is_Symbol:
            used_symbols.append(term)
        else:
            used_symbols.extend(retrieve_used_symbols(term))
    return used_symbols

当你有混合符号时,后者会派上用场:

>>> crazy_expr = expr + 10*y-2*z
>>> crazy_expr
a[1] + 4*a[4] + 10*y - 2*z + 5
>>> used_symbols = retrieve_used_symbols(crazy_expr)
>>> replacements = [(n, x**(enu+1)) for enu,n in enumerate(used_symbols)]
>>> Poly(crazy_expr.subs(replacements)).all_coeffs()
[4, -2, 1, 10, 5]
>>> list(reversed(used_symbols))
[a[4], z, a[1], y]

对于IndexedBase对象,更简单:

coeffs = [expr.coeff(x[i]) for i in range(10)]

但是您仍然需要添加常数项,如您所说,您可以通过通配符替换获得该常数项:

ind = Wild('i')
constant_term = expr.replace(x[ind], 0)

{根据@(Oliver W.) 的要求}

给出

>>> x = IndexedBase('x')
>>> eqs = 4*x[0] + 5*x[5] + 1, x[1] - x[2]
>>> v = list(ordered(Tuple(*eqs).atoms(Indexed)))

可以这样做

>>> [[eq.coeff(vi) for vi in v] + [eq.as_coeff_Add()[0]] for eq in eqs] 
[[4, 0, 0, 5, 1], [0, 1, -1, 0, 0]]

但是其中大部分可以通过矩阵方法获得jacobian。但是要使用它你必须用符号替换 x[i] (因为 diff 只适用于函数是符号,IIRC):

>>> d = [Dummy() for vi in v]
>>> z = dict(zip(d, [0]*len(d)))
>>> m = Matrix([eq.xreplace(dict(zip(v, d))) for eq in eqs])
>>> m.jacobian(d)
Matrix([
[4, 0,  0, 5],
[0, 1, -1, 0]])
>>> m.subs(z)
Matrix([
[1],
[0]])