高效地评估 80 个变量的多项式，其中包含 80,000 个项
Efficiently evaluate polynomial of 80 variables with 80,000 terms
我有一个包含 80 个特征（列）和 4000 个实例（行）的 CSV 文件。我想将此文件转换为更高维度（内核 3）。新文件应该有 88560 列和 4000 行。我正在尝试做一个只有 3 个值的试验。我首先准备了带有符号 a、b、c 的 header，然后尝试输入值但无法用值替换字符串。代码如下。
import csv
fw = open('polyexp.csv', "w")
a,b,c= sympy.symbols("a b c")
formula = (a+b+c) ** 3
poly = str(formula.expand())
ls = poly.split('+')
print >> fw, ls
value1 = ls[0]
value2 = ls[1]
value3 = ls[3]
a= 1
b=2
c =3
print value1,value2,value3
这给出输出 a**3 3*a**2*b 3*a*b**2 而不是 1 6 12。
请问还有其他方向吗！
我正在发布我修改后的代码，它正在运行，但是有什么优雅的方法可以做到这一点，特别是部分 1) 将 80 个特征转换为符号 2) 替换它们现在硬编码的值
import csv, sympy
import sympy as sym
fw = open('polyexp.csv', "w")
fr = open('polyval.csv', "r")
a,b,c, d, e= sympy.symbols("a b c d e")
formula = (a+b+c+d+e) ** 3
poly = str(formula.expand())
ls = poly.split('+')
for i in range(len(ls)):
ls_val[i] = sym.sympify(ls[i])
for line in fr:
tup = line.split(",")
csvrow = list()
for key in range(len(ls_val)):
new = ls_val[key].subs([(a, tup[0]),(b,tup[1]),(c,tup[2]),(d,tup[3]),(e,tup[4])])
csvrow.append(new)
print >> fw, csvrow
在您的代码中，value1,value2,value3 只是字符串，因为您使用了 str(formula.expand())。你需要再次从他们那里做出 sympy 表达：
value1 = sympify(ls[0])
value2 = sympify(ls[1])
value3 = sympify(ls[3])
然后，您可以将其中的符号替换为您的实际值以获得数值结果：
print value1.subs([(a, 1), (b, 2), (c, 3)])
print value2.subs([(a, 1), (b, 2), (c, 3)])
print value3.subs([(a, 1), (b, 2), (c, 3)])
不要像这样将值 1, 2, 3 分配给符号 a, b, c：a = 1。这样你就失去了符号并且 a 变成了一个整数 1 而不是具有该值的符号。
另见 tutorial。
这是一个生成 80 个符号的版本，但仍然相当慢。
根据文档，它使用 sympy.lambdify() 比 subs()/evalf() 快 ~50 倍。
import sympy as sym
fw = open('polyexp.csv', "w")
fr = open('polyval.csv', "r")
flabel = open('polylabel.csv', "r")
N = 80
symbolnames = ['a'+str(i) for i in range(N)]
symbols = [sym.Symbol(name) for name in symbolnames] # Generate symbols 'a0', 'a1', ...
formula = sym.Add(*symbols) ** 3
poly = str(formula.expand())
terms = [sym.sympify(term) for term in poly.split('+')]
funcs = [sym.lambdify(symbolnames, term, 'math') for term in terms]
for line in fr:
label = flabel.readline().rstrip()
values = [float(s) for s in line.split(",")]
csvrow = [func(*values) for func in funcs]
print >> fw, label + ',' + ','.join(map(str, csvrow)
最后，这是一个使用 numpy 的非常快速的版本，其中 lambda 函数的一次调用计算输出文件的一整列。由于这种按列处理，整个结果数组在写出之前必须保存在内存中。如果你没有足够的内存，你可以把结果列写成行，给你一个转置的输出文件。
import numpy as np
import sympy as sym
fw = open('polyexp.csv', "w")
flabel = open('polylabel.csv', "r")
N = 80
symbolnames = ['a{:02}'.format(i) for i in range(N)]
# Try to read polynomial expansion from cache file (saves time).
polycachefile = 'poly{}.txt'.format(N)
poly = ''
try:
with open(polycachefile) as f:
poly = f.readline().strip()
except IOError:
poly = ''
if poly.count('+') > 0:
# Use cached polynomial expansion.
else:
# Calculate and save polynomial expansion.
symbols = [sym.Symbol(name) for name in symbolnames]
formula = sym.Add(*symbols) ** 3
poly = str(formula.expand())
with open(polycachefile, 'w') as f:
f.write(poly)
terms = poly.split('+')
# Read input file into val.
val = np.genfromtxt('polyval.csv', delimiter=',', \
autostrip=True, dtype=np.float32) # <-- adjust data type (i.e. precision) for calculations here!
nrows, ncols = val.shape
assert ncols == N
# Prepare efficient access to columns of val.
columns = val.T[:, :]
colsbyname = dict(zip(symbolnames, columns))
symbolnamesset = set(symbolnames)
# Create result array.
exp = np.zeros([nrows, len(terms)], dtype=val.dtype)
# Calculate result, column by column.
for i, term in enumerate(terms):
term = term.strip()
subterms = set(term.split('*'))
usedsyms = subterms & symbolnamesset
func = sym.lambdify(usedsyms, term, 'math')
exp[:, i] = func(*(colsbyname[s] for s in usedsyms))
# Write result file.
rowfmt = ','.join(['%.7g'] * len(terms)) # <-- adjust output precision here!
for row in exp:
label = flabel.readline().rstrip()
print >> fw, label + ',' + rowfmt % tuple(row)
fw.close()
性能：在我的 Core i3 上，计算需要 35 秒，写入结果文件需要 8 分钟才能达到 17 位精度（或 3 分钟才能达到 7 位精度）。
我有一个包含 80 个特征（列）和 4000 个实例（行）的 CSV 文件。我想将此文件转换为更高维度（内核 3）。新文件应该有 88560 列和 4000 行。我正在尝试做一个只有 3 个值的试验。我首先准备了带有符号 a、b、c 的 header，然后尝试输入值但无法用值替换字符串。代码如下。
import csv
fw = open('polyexp.csv', "w")
a,b,c= sympy.symbols("a b c")
formula = (a+b+c) ** 3
poly = str(formula.expand())
ls = poly.split('+')
print >> fw, ls
value1 = ls[0]
value2 = ls[1]
value3 = ls[3]
a= 1
b=2
c =3
print value1,value2,value3
这给出输出 a**3 3*a**2*b 3*a*b**2 而不是 1 6 12。
请问还有其他方向吗！
我正在发布我修改后的代码，它正在运行，但是有什么优雅的方法可以做到这一点，特别是部分 1) 将 80 个特征转换为符号 2) 替换它们现在硬编码的值
import csv, sympy
import sympy as sym
fw = open('polyexp.csv', "w")
fr = open('polyval.csv', "r")
a,b,c, d, e= sympy.symbols("a b c d e")
formula = (a+b+c+d+e) ** 3
poly = str(formula.expand())
ls = poly.split('+')
for i in range(len(ls)):
ls_val[i] = sym.sympify(ls[i])
for line in fr:
tup = line.split(",")
csvrow = list()
for key in range(len(ls_val)):
new = ls_val[key].subs([(a, tup[0]),(b,tup[1]),(c,tup[2]),(d,tup[3]),(e,tup[4])])
csvrow.append(new)
print >> fw, csvrow
在您的代码中，value1,value2,value3 只是字符串，因为您使用了 str(formula.expand())。你需要再次从他们那里做出 sympy 表达：
value1 = sympify(ls[0])
value2 = sympify(ls[1])
value3 = sympify(ls[3])
然后，您可以将其中的符号替换为您的实际值以获得数值结果：
print value1.subs([(a, 1), (b, 2), (c, 3)])
print value2.subs([(a, 1), (b, 2), (c, 3)])
print value3.subs([(a, 1), (b, 2), (c, 3)])
不要像这样将值 1, 2, 3 分配给符号 a, b, c：a = 1。这样你就失去了符号并且 a 变成了一个整数 1 而不是具有该值的符号。
另见 tutorial。
这是一个生成 80 个符号的版本，但仍然相当慢。
根据文档，它使用 sympy.lambdify() 比 subs()/evalf() 快 ~50 倍。
import sympy as sym
fw = open('polyexp.csv', "w")
fr = open('polyval.csv', "r")
flabel = open('polylabel.csv', "r")
N = 80
symbolnames = ['a'+str(i) for i in range(N)]
symbols = [sym.Symbol(name) for name in symbolnames] # Generate symbols 'a0', 'a1', ...
formula = sym.Add(*symbols) ** 3
poly = str(formula.expand())
terms = [sym.sympify(term) for term in poly.split('+')]
funcs = [sym.lambdify(symbolnames, term, 'math') for term in terms]
for line in fr:
label = flabel.readline().rstrip()
values = [float(s) for s in line.split(",")]
csvrow = [func(*values) for func in funcs]
print >> fw, label + ',' + ','.join(map(str, csvrow)
最后，这是一个使用 numpy 的非常快速的版本，其中 lambda 函数的一次调用计算输出文件的一整列。由于这种按列处理，整个结果数组在写出之前必须保存在内存中。如果你没有足够的内存，你可以把结果列写成行，给你一个转置的输出文件。
import numpy as np
import sympy as sym
fw = open('polyexp.csv', "w")
flabel = open('polylabel.csv', "r")
N = 80
symbolnames = ['a{:02}'.format(i) for i in range(N)]
# Try to read polynomial expansion from cache file (saves time).
polycachefile = 'poly{}.txt'.format(N)
poly = ''
try:
with open(polycachefile) as f:
poly = f.readline().strip()
except IOError:
poly = ''
if poly.count('+') > 0:
# Use cached polynomial expansion.
else:
# Calculate and save polynomial expansion.
symbols = [sym.Symbol(name) for name in symbolnames]
formula = sym.Add(*symbols) ** 3
poly = str(formula.expand())
with open(polycachefile, 'w') as f:
f.write(poly)
terms = poly.split('+')
# Read input file into val.
val = np.genfromtxt('polyval.csv', delimiter=',', \
autostrip=True, dtype=np.float32) # <-- adjust data type (i.e. precision) for calculations here!
nrows, ncols = val.shape
assert ncols == N
# Prepare efficient access to columns of val.
columns = val.T[:, :]
colsbyname = dict(zip(symbolnames, columns))
symbolnamesset = set(symbolnames)
# Create result array.
exp = np.zeros([nrows, len(terms)], dtype=val.dtype)
# Calculate result, column by column.
for i, term in enumerate(terms):
term = term.strip()
subterms = set(term.split('*'))
usedsyms = subterms & symbolnamesset
func = sym.lambdify(usedsyms, term, 'math')
exp[:, i] = func(*(colsbyname[s] for s in usedsyms))
# Write result file.
rowfmt = ','.join(['%.7g'] * len(terms)) # <-- adjust output precision here!
for row in exp:
label = flabel.readline().rstrip()
print >> fw, label + ',' + rowfmt % tuple(row)
fw.close()
性能：在我的 Core i3 上，计算需要 35 秒，写入结果文件需要 8 分钟才能达到 17 位精度（或 3 分钟才能达到 7 位精度）。