打印 syms / matlabFunction 慢答案

【问题标题】：Printing syms / matlabFunction slow打印 syms / matlabFunction 慢
【发布时间】：2015-10-16 16:06:31
【问题描述】：

我在尝试使符号替换更快时遇到了很多麻烦 - 也就是说，用符号表达式替换变量并得到一个双精度。

我正在创建一个复杂的函数 f，并计算它的雅可比 df。这以合理的速度进行，我可以将其保存到文件中。但是当我尝试使用 matlabFunction 甚至 disp 或 fprintf 时，系统会挂起并且无法继续进行（即使 matlabFunction 设置为未优化）。这是一个主要问题，因为我需要能够进行相当快的替换。

f 向量是 24 个元素，雅可比是 24 x 78（不过这里只显示了 70 个变量，所以可以压缩到 70 列；不过我怀疑这是问题所在）。

我还知道 f 和 df 的某些元素在单独访问时很简单并且工作正常，但是 f 和 df 的某些更复杂的元素无法显示。我想它们很长，但由于它们计算得很好，所以对我来说它们不能转换为 matlabFunction 或显示是没有意义的。

更奇怪的是，我可以用 in 替换我所有的符号变量，但是完全替换的 f 向量的最终显示（例如，通过 disp）或转换为 double（通过 double()）似乎需要永远.

如果你想玩弄 .mat 文件，你可以得到它here（filedropper 链接，它是 288kb）。如何在合理的时间内写出这个文件？

【问题讨论】：

我正在写一个格式化的评论作为答案，尝试simplify你的函数，并通过让 matlab 吃掉所有的内存（失去我的答案）迅速终止了我的会话。我怀疑这可能是在你的挂起。无论如何：你的f(13) 太可怕了，超过 25k 个字符长（这是 matlab 的最大值），字面常量大约是 1e33 和像 BASE_ORIGIN_Y 这样的东西。就不能以某种方式避开这个怪物吗？我认为在目前的形式下，绝对没有办法将其转换为易于处理的函数，您可以将单个值替换为它本身就是一个奇迹。
它是矩阵乘法和加法的结果。我一直尝试使用简化命令对其进行简化，但这非常慢。我曾经成功地将 f 写入一个未优化的文件（df 没有运气）——花了 2 个小时——但随后评估它花了 0.8 秒，这太慢了。我需要能够在大约 0.02 秒内执行评估。我现在没有看到用先验知识简化它的分析方法，所以我正在寻找关于使 syms 更快地工作的方法的建议，或者可以完成工作的不同工具（框）。
对matlab来说是个无望的案例。切换到 Maple 或 Mathematica 用于 CAS 目的，甚至尝试 SymPy
您认为 Sympy 或 Mathematica 能够处理这个问题吗？为什么你这么想？你能提供一些实验或其他东西来验证吗？如果该方法可行，我可以将我的数据在那里洗牌，然后再洗牌。

标签： matlab symbolic-math mupad

【解决方案1】：

注意：从您的评论的角度来看，我专注于这个问题：

它是矩阵乘法和加法的结果。我一直尝试使用简化命令对其进行简化，但这非常慢。我曾经成功地将 f 写入一个未优化的文件（df 没有运气）——花了 2 个小时——但随后评估它花了 0.8 秒，这太慢了。我需要能够在大约 0.02 秒内执行评估。

我开始查看您的f 中的元素，直到f(12) 都很简单。然而，f(13) 释放了地狱：

>> inp.f(13)

ans =

(2289*l4)/100 - (11371197146449238679*l3)/8112963841460668169578900514406400 - (2289*l2)/100 + (11371197146449238679*l5)/8112963841460668169578900514406400 - (2289*l8)/100 - (11371197146449238679*l9)/8112963841460668169578900514406400 + (2289*l10)/100 + (11371197146449238679*l11)/8112963841460668169578900514406400 - (2289*l14)/100 - (11371197146449238679*l15)/8112963841460668169578900514406400 + (2289*l16)/100 + (11371197146449238679*l17)/8112963841460668169578900514406400 - (2289*l20)/100 - (11371197146449238679*l21)/8112963841460668169578900514406400 + (2289*l22)/100 + (11371197146449238679*l23)/8112963841460668169578900514406400 - (2289*l26)/100 - (11371197146449238679*l27)/8112963841460668169578900514406400 + (2289*l28)/100 + (11371197146449238679*l29)/8112963841460668169578900514406400 - (2289*l32)/100 - (11371197146449238679*l33)/8112963841460668169578900514406400 + (2289*l34)/100 + (11371197146449238679*l35)/8112963841460668169578900514406400 - h1*(((cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))*(cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2)) + (sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2))*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) - cos(x5/2)^2*cos(x6/2)*sin(x6/2))*(((x17*(cos(x4/2)*cos(x5/2)*(cos(x6/2)*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) + sin(x6/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2))) - cos(x5/2)*sin(x4/2)*(cos(x6/2)*(cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2)) + sin(x6/2)*(cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2)))))/2 - (x18*(cos(x4/2)^2*cos(x5/2)^2 + cos(x5/2)^2*sin(x4/2)^2 + sin(x5/2)^2))/2 + (x16*(sin(x5/2)*(cos(x5/2)^2*cos(x6/2)^2 + cos(x5/2)^2*sin(x6/2)^2 + sin(x5/2)^2) + cos(x5/2)*sin(x4/2)*(cos(x5/2)*sin(x4/2)*sin(x5/2) + cos(x5/2)*cos(x6/2)*(cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2)) - cos(x5/2)*sin(x6/2)*(cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))) + cos(x4/2)*cos(x5/2)*(cos(x4/2)*cos(x5/2)*sin(x5/2) - cos(x5/2)*cos(x6/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + cos(x5/2)*sin(x6/2)*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)))))/2 - (x19*cos(x5/2)*sin(x4/2))/2)*((LEG_MASS*((cos(x7/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + cos(x5/2)*cos(x6/2)*sin(x7/2))*((cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))*(x2/2 - BASE_ORIGIN_Z*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) - (cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*sin(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4)) - (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*(x1/2 + BASE_ORIGIN_Z*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*cos(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4)) + cos(x5/2)*sin(x4/2)*(x3/2 - sin(x5/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4) + BASE_ORIGIN_Z*cos(x4/2)*cos(x5/2) - cos(x5/2)*sin(x4/2)*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2))) - (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*((sin(x5/2)*sin(x7/2) - cos(x4/2)*cos(x5/2)*cos(x7/2))*(x3/2 - sin(x5/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4) + BASE_ORIGIN_Z*cos(x4/2)*cos(x5/2) - cos(x5/2)*sin(x4/2)*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2)) - (cos(x7/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + cos(x5/2)*cos(x6/2)*sin(x7/2))*(x1/2 + BASE_ORIGIN_Z*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*cos(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4)) + (cos(x7/2)*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) - cos(x5/2)*sin(x6/2)*sin(x7/2))*(x2/2 - BASE_ORIGIN_Z*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) - (cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*sin(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4))))*(sin(x7/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) - cos(x5/2)*cos(x6/2)*cos(x7/2)) - LEG_MASS*((sin(x7/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) - cos(x5/2)*cos(x6/2)*cos(x7/2))*((cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))*(x2/2 - BASE_ORIGIN_Z*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) - (cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*sin(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4)) - (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*(x1/2 + BASE_ORIGIN_Z*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*cos(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4)) + cos(x5/2)*sin(x4/2)*(x3/2 - sin(x5/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4) + BASE_ORIGIN_Z*cos(x4/2)*cos(x5/2) - cos(x5/2)*sin(x4/2)*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2))) + (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*((cos(x7/2)*sin(x5/2) + cos(x4/2)*cos(x5/2)*sin(x7/2))*(x3/2 - sin(x5/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4) + BASE_ORIGIN_Z*cos(x4/2)*cos(x5/2) - cos(x5/2)*sin(x4/2)*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2)) + (sin(x7/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) - cos(x5/2)*cos(x6/2)*cos(x7/2))*(x1/2 + BASE_ORIGIN_Z*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*cos(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4)) - (sin(x7/2)*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) + cos(x5/2)*cos(x7/2)*sin(x6/2))*(x2/2 - BASE_ORIGIN_Z*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) - (cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*sin(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4))))*(cos(x7/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + cos(x5/2)*cos(x6/2)*sin(x7/2)) + LEG_MASS*(cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*((cos(x7/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + cos(x5/2)*cos(x6/2)*sin(x7/2))*((cos(x7/2)*sin(x5/2) + cos(x4/2)*cos(x5/2)*sin(x7/2))*(x3/2 - sin(x5/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4) + BASE_ORIGIN_Z*cos(x4/2)*cos(x5/2) - cos(x5/2)*sin(x4/2)*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2)) + (sin(x7/2)*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) - cos(x5/2)*cos(x6/2)*cos(x7/2))*(x1/2 + BASE_ORIGIN_Z*(sin(x4/2)*sin(x6/2) + cos(x4/2)*cos(x6/2)*sin(x5/2)) + (cos(x4/2)*sin(x6/2) - cos(x6/2)*sin(x4/2)*sin(x5/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + cos(x5/2)*cos(x6/2)*(BASE_ORIGIN_X - BASE_LINK_EXTENTS_X/4)) - (sin(x7/2)*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) + cos(x5/2)*cos(x7/2)*sin(x6/2))*(x2/2 - BASE_ORIGIN_Z*(cos(x6/2)*sin(x4/2) - cos(x4/2)*sin(x5/2)*sin(x6/2)) - (cos(x4/2)*cos(x6/2) + sin(x4/2)*sin(x5/2)*sin(x6/2))*(BASE_LINK_EXTENTS_Y/2 - BASE_ORIGIN_Y + LEG_LINK_EXTENTS_Y/2) + [...] ... Output truncated.  Text exceeds maximum line length of 25,000 characters for Command Window display.

（实际输出被截断以适应 SO 的 30k 字符限制，但你明白了）

我敢打赌matlabFunction 的解析器不适合这种量级的输入。里面还有一些奇怪的东西：比如8e33顺序的整数字符串文字。

所以我仔细研究了你的功能。幸运的是，您可以将函数转换为字符串，然后处理这些字符串，这只会占用 CPU 时间而不占用内存。

预处理：

for k=1:24
   fstring2{k}=char(inp.f(k));
end

函数长度：

>> cellfun(@length,fstring2)

ans =

  Columns 1 through 12

          11          11          11          11          11          11          11          11          11          11          11          11

  Columns 13 through 24

     2301006     2300241     2299996     8425640     8416273     8424306     1375443     1305245     1302440     1237876     1381084     1310884

休斯顿，我们遇到了问题。

这些符号函数的巨大野兽破坏了matlabFunction的解析器，或者更有可能的是，您在操作过程中耗尽了内存。当我尝试 simplify f(13) 时，我确实做到了，在几秒钟内丢失了 8 GB 的一半。

作为概念验证，我试图模拟您的函数所涉及的计算工作。我检查了f(13)（第一只野兽）。有关所涉及操作的一些信息：

>> length(strfind(char(inp.f(13)),'*'))

ans =

      134710

>> length(strfind(char(inp.f(13)),'+'))

ans =

       36932

>> length(strfind(char(inp.f(13)),'-'))

ans =

       26855

>> length(strfind(char(inp.f(13)),'/'))

ans =

      183380

>> length(strfind(char(inp.f(13)),'ln'))

ans =

     0

>> length(strfind(char(inp.f(13)),'exp'))

ans =

     0

>> length(strfind(char(inp.f(13)),'cos'))

ans =

       78700

>> length(strfind(char(inp.f(13)),'sin'))

ans =

       84142

我尝试对涉及相似数量操作的模拟计算进行计时：

x=zeros(36000,1);
tic;
for k=1:36000
   x(k)=(((sin(sin(((cos(cos(3.1+2.1)*3.1)*6.1)*5.1)*9.1)/4.1)/3.1)/6.1)/5.1)/8.1;
end
toc;

Elapsed time is 0.010895 seconds.

这涉及36000 加法、144000 乘法、180000 除法和72000 对sin 和cos 的调用。

现在，如果我们假设这是一个正确的大致数字，并且如果我们假设您的函数具有类似的操作分布，那么您正在查看函数的 40080434 字符，相当于 17 个 f(13)单位。这表明即使您可以转换为适当的 matlab 函数，您的运行时仅调用 f（我们根本没有查看 df）至少需要 0.1-0.2秒。

由于您的问题的性质，我不确定是否有解决方法。我可能会尝试在python 中使用sympy 做同样的事情，您还可以转换为lambda（python 等效于匿名函数）用于数值计算。如果 that 会成功，那么至少你可以尽快使用你的函数。

更新

在发布了我不太乐观的答案后，我相信我已经成功地将您的函数转换为匿名函数。它很脏，但它似乎工作。

首先将函数转换为上述字符串，然后使用symvar 提取变量名。然后使用这些函数名创建一个函数定义；不幸的是，我只能使用eval 破解它。应该有一种更优雅的方式，但无论如何我们对可实现的运行时感兴趣。

varcell=symvar(fstring2{13}); %variables of inp.f(13)
vars2=strcat(varcell,','); %add a comma to each var
vars3=[vars2{:}]; %put them into a single string
vars3=vars3(1:end-1); %remove trailing comma

f13=eval(['@(' v3 ') ' fstring2{13}]); %this is your numeric function

转换很麻烦，但匿名函数的实际构造很快，而且不会太占用内存。虚拟运行时：

>> tic; ftry(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58), toc

ans =

   1.1417e+06

Elapsed time is 0.069252 seconds.

它可以变得更加用户友好，例如允许在函数中进行数组操作，或者将所有 58 个输入作为单个数组输入传递。但是你的运行时间是一样的。这只是一个功能，您大约有 17 个。您可能永远无法获得您希望的加速。

（无论如何，我确实开始得到

线程“AWT-EventQueue-0”java.lang.OutOfMemoryError 中的异常：Java 堆空间

整个磨难之后的错误，所以它的成功也可能取决于你对“成功”的定义；）

【讨论】：

我能够了解您使用 matlabFunction 的情况，但它仅适用于 f。这可能为我提供了一种 df 的交互方式，值得赞赏。这仍然不是我需要的速度，所以如果有一种方法可以在这里运行简化，那就太好了。但这似乎太贵了，而且没有办法加快速度。为了赏金，我将为您在这里的所有辛勤工作给您 50 分。但我不会将其标记为已选中。我将继续努力加快速度并尽快报告。
@user650261 非常感谢，您说得有道理:) 您可以尝试在 python 中使用 sympy 与 numpy 做同样的事情吗？
@user650261 仅供参考：我将f(13) 作为字符串，将其打印到文件中并用它来尝试 python。如果我定义一个基于 numpy 的数字 lambda（python 等效的匿名函数），或者如果我将其定义为符号函数，然后让 sympy 将其转换为 lambda：这两种情况都需要 0.46 秒来计算我的上面的最后一个示例，无论数字/符号来源如何。因此，即使从 sympy 转到 numpy 也可以，但它似乎 太慢了...
Daek 感谢您对此进行调查。该项目基于 MATLAB，因此我无法将所有代码都转换为 Python。如果我能以某种方式获得一些性能提升，那么我换用 Python 才值得。当我有机会时，我想尝试的一件事是看看 sympy 的简化方法是否可以在比 MATLAB 更合理的时间内工作，并看看在那里提供了哪些性能加速。
@user650261 我也开始研究这个问题。在一个小时和 5GB 的内存中，f(13) 没有得到任何结果，所以我不知道如果它最终能找到它的效率会有多高。