Login/Register
Printing syms / matlabFunction slow
Asked Answered
matlabsymbolic-mathmupad
I

1

5

I am having a lot of trouble trying to make symbolic substitution go faster - that is, substituting in for variables in a symbolic expression and getting out a double.

I am creating a complicated function f, and calculating its jacobian df. This goes at a reasonable pace, and I can save it to a file just fine. But when I try to use matlabFunction or even disp or fprintf, the system hangs and is unable to proceed further (even when matlabFunction is set to unoptimized). This is a major problem as I need to be able to do reasonably fast substitution.

The f vector is 24 elements, and the Jacobian is 24 x 78 (there are only 70 variables that show up here though, so this could be compressed down to 70 columns; I doubt this is the issue though).

I also know that there are certain elements of f and df which are simple and work fine when accessed individually, but certain, more complicated elements of f and df cannot be displayed. I imagine they are very long, but since they are calculated just fine, it doesn't make sense to me that they cannot be transformed into a matlabFunction or displayed.

Even more bizarrely, I can substitute in for all my symbolic variables, but then the final display of the fully substituted f vector (say, via disp), or the conversion to double (through double( )) seems to take forever.

If you want to play around with the .mat file, you can get it here (filedropper link, it's 288kb). What can I do to write out this file in a reasonable amount of time?

Iggie answered 16/10, 2015 at 16:6 Comment(4)
I was writing a formatted comment as answer, tried to simplify your function and promptly killed my session by letting matlab eat all the memory (losing my answer). I suspect this might be behind your hangs. Anyway: your f(13) is horrifying, more than 25k characters long (which is the max for matlab) with literal constants of around 1e33 and things like BASE_ORIGIN_Y in there. Can't this monster be avoided in some way? I think in its current form there is absolutely no way to convert it to a tractable function, it's a miracle in itself that you can substitute individual values into it. – Dowlen
It is the result of a matrix multiplication and addition. I have tried simplifying it along the way, using the simplify command, but that's been incredibly slow. I have succeeded once in writing f to a file unoptimized (no luck for df) - it took 2 hours - but then evaluating it took 0.8 seconds, which is too slow. I need to be able to perform evaluation in about 0.02 seconds. I don't see an analytical way right now of simplifying it with a priori knowledge, so I'm looking for suggestions of either ways to make syms to work faster, or a different tool(box) that can get the job done. – Iggie
It is an hopeless case for matlab. Switch to Maple or Mathematica for CAS purposes or even try SymPy – Paste
Do you think Sympy or Mathematica will be able to handle this? Why do you think that? Can you provide some experiment or something to validate that? If that method will work, I can shuffle my data there and then shuffle it back. – Iggie
V
7

Note: I concentrated on the problem from the viewpoint of your comment:

It is the result of a matrix multiplication and addition. I have tried simplifying it along the way, using the simplify command, but that's been incredibly slow. I have succeeded once in writing f to a file unoptimized (no luck for df) - it took 2 hours - but then evaluating it took 0.8 seconds, which is too slow. I need to be able to perform evaluation in about 0.02 seconds.

I started looking at the elements in your f, and it was simple up to f(12). However, f(13) unleashed hell:

>> 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.

(actual output truncated to fit in 30k character limit of SO, but you get the deal)

I'd wager the parser of matlabFunction was not meant for inputs of this magnitude. There are also some weird things in there: like integer string literals on the order of 8e33.

So I took a closer look at your function. Fortunately you can convert your functions to strings, and work on those, which is only heavy on CPU time but not on memory.

Preproc:

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

Function lengths:

>> 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

Houston, we've had a problem.

These massive beasts of symbolic functions break the parser of matlabFunction, or what's more probable, you run out of memory during the operation. I sure did when I tried to simplify f(13), lost the better half of 8 GB within seconds.

Just as a proof-of-concept I tried to mock the computational effort involved in your functions. I inspected f(13) (the first beast). Some info about the operations involved:

>> 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

I tried to time a mock computation involving a similar number of operations:

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.

This involves 36000 additions, 144000 multiplications, 180000 divisions and 72000 calls to sin and cos, each.

Now, if we assume that this is a correct ballpark figure, and if we assume that your functions have a similar distribution of operations, then you're looking at 40080434 characters of functions, which is 17 equivalent f(13) units. This suggests that even if you could convert to a proper matlab function, your runtime to just calling f (and we haven't looked at df at all) would take at least 0.1-0.2 seconds.

Due to the nature of your problem, I'm not sure there's a way around it. I would probably try doing the same using sympy in python, there you can also convert to a lambda (the python equivalent of an anonymous function) for use in numerical computations. If that would succeed, then at least you could use your functions as quickly as possible.

UPDATE

After posting my less optimistic answer, I believe I've succeeded in converting your function to an anonymous one. It's dirty, but it seems to work.

First you convert your function to a string as above, then use symvar to extract the variable names. Then you create a function definition using these function names; unfortunately I could only hack it up using eval. There should be a more elegant way, but anyway we're interested in the achievable runtimes.

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

The conversion is nasty, but the actual construction of the anonymous function is quick and not too hard on memory. Dummy runtime:

>> 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.

It could be made much more user-friendly, for instance by allowing array operations in the function, or passing all 58 inputs as a single array input. But yor runtime will be the same. And this is just one function, and you have roughly 17 of these. You might never get the speed-up you're hoping for.

(And anyway, I did start getting

Exception in thread "AWT-EventQueue-0" java.lang.OutOfMemoryError: Java heap space

errors after this whole ordeal, so its success might also depend on your definition of "success";)

Various answered 20/10, 2015 at 11:46 Comment(5)
I was able to get about as far as you were using matlabFunction, but it only worked for f. This probably gives me an interactive way for df, which is appreciated. It's still not the speed I need, so it would be nice if there were a way to run a simplify here. But it seems like that's too expensive and there's no way to speed this up. For the sake of the bounty, I'm going to give you the 50 points for all your hard work here. But I'm not going to mark this as checked. I will continue trying to make this faster and report back soon. – Iggie
@Iggie thank you very much, you have a fair point:) Is it an option for you to try and do the same in python using sympy vs numpy? – Dowlen
@Iggie just FYI: I took f(13) as a string, printed it into a file and used that to try python. If I define a numpy-based numeric lambda (the python equivalent of an anonymous function), or if I define it as a symbolic function and then let sympy convert it to a lambda: both cases take 0.46 s to compute my last example above, regardless of the numeric/symbolic origin. So even though it would be OK to go from sympy to numpy, it just seems soooo slow... – Dowlen
Daek Thanks for looking into that. This project is based in MATLAB, so I can't convert all my code to Python. It would only be worth me switching to Python if I could get some performance boost in some way. The one thing I want to try when I get a chance is to see if sympy's simplify methods might work in a more reasonable amount of time than MATLAB's, and see what performance speedup that gives there. – Iggie
@Iggie I started looking into that as well. In an hour and 5GB of memory it didn't get anywhere with f(13), so I don't know how efficient the result would be if it would find it eventually. – Dowlen

© 2022 - 2024 — McMap. All rights reserved.