how to make symbolic computations fast?

Question

hosein Javan 2020-8-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/580290-how-to-make-symbolic-computations-fast

编辑： Pavel Holoborodko 2020-8-29

hello, I'm working on some calculation that requires great precision at very large or small values.

as I have searched the only way is symbolic "vpa". this way is more precise with a maximum of 2^29 digits but is really slow. for example consider that we are going to estimate bessel function of first kind.

% use double precision data (not precise but fast)
tic
besseli(1,100)
toc
>>ans = 1.068369390338163e+42
Elapsed time is 0.000295 seconds.
% use variable precision arithmetic data (precise but not fast)
tic
vpa(besseli(1,sym(100)),42)
toc
>>ans = 1068369390338162481206145763224295265446122.0
Elapsed time is 0.701682 seconds.

But I was wondering if there is a way somehow to store data and use interpolation function to speed it up. I list my questions as following:

is there any way to speed up symbolic computation?
is there memory-based method that reduces cpu consumption like "look-up tables and pre strored data interpolation" for symbolic numbers?
if this way is some how inefficient, what other best method is proposed?

4 个评论
显示 2更早的评论隐藏 2更早的评论

Steven Lord 2020-8-18

What are you doing that you need (or think you need) such high precision? What problem are you trying to solve? Maybe there's a way to solve that problem that doesn't involve increased precision.

hosein Javan 2020-8-18

Steven Lord. please read the comments. thanks.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2020-8-16

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/580290-how-to-make-symbolic-computations-fast#answer_480756

编辑：John D'Errico 2020-8-16

在 MATLAB Online 中打开

No. Your conclusion the only way is VPA is simply wrong. You can use my HPF toolbox instead. It is available for free download from the File Exchange. However, I'm sorry, but interpolation simply makes virtually no sense in any such context.

https://www.mathworks.com/matlabcentral/fileexchange/36534-hpf-a-big-decimal-class

Using both tools to compute 1000 digits for e = exp(1) using a simple approximation

https://en.wikipedia.org/wiki/List_of_representations_of_e

I might do this:

DefaultNumberOfDigits 1000
tic
eappr = hpf(2.5); % Be careful. I know 2.5 is exactly representable as a double.
den = hpf(2);
for k = 1:150
  den = den*(3*k)*(3*k+1)*(3*k+2);
  eappr = eappr + hpf((3*k + 2)^2 +1)/den;
end
toc
exp(hpf(1)) - eappr
ans.Exponent
Elapsed time is 0.693297 seconds.
ans =
6.e-1009

Now for vpa...

digits 1000
tic
eappr = vpa(2.5);
den = vpa(2);
for k=1:270
  den = den*(3*k)*(3*k+1)*(3*k+2);
  eappr = eappr + sym((3*k + 2)^2 +1)/den;
end
toc
exp(vpa(1)) - eappr
Elapsed time is 0.516161 seconds.
ans =
0.0

150 terms from that series are sufficient to compute 1000 digits for exp(1).

In both cases, they are correct to the full precision required. VPA beat HPF in time by a little, but not a lot.

I will add that I used HPF to compute a million digits of pi once.

I would also add that much of the time, people are incorrrect in thinking they need to use high precision arithmetic to do computations. It is often the case they could have used logs to do essentially everything they really needed. But that is not always the case.

14 个评论
显示 12更早的评论隐藏 12更早的评论

hosein Javan 2020-8-16

hello, Mr. D'Errico, it's nice to speak with u. I have heard of ur work from Mr. Roberson. your work is great. but I assume it is not suitable for my problem. But I am dealing with a problem that has taken a lot of time of me about one month just because of software problems. well in order to shorten the topic I simply couldn't tell everything.

this program, I'm writing needs a speed way quicker than this, due to Large amount of data in multiple loops for future optimization and design purposes. moreover I have functions that are really complicated (like integral of product of bessel functions), that functions requires 8 seconds to just calculate for x=1.

there's also a need for computing very low values like 1e-400, does ur toolbox HPF and VPI support that too?

I've also tried methods for avoiding large computations, like scaling, unit changing but I found that scaling is not even an idea because for example "1e400/1e400=inf~=1", and unit changing does not work for me because there's a factor that when I decrease the unit, it increases and vice versa.

I've also found that there's a software developer "advanpix" that has managed to write a code that supports multi-precision values like quadruple but that is not for free.

-------------------------------------------

if you have time, I would like to explain a brief explanation of what I'm about to do. I'm trying to solve a first order bessel equation which is part of a more complicated PDE problem, whose solution is as following form:

now we have a bunch of these functions and not one, so there are for example 10 integral coefficients. now when I apply the boundary conditions, a set of linear equations are made by which I can calculate c1, c2,... . in the form "A*C=B".

now as x is near zero, I1(x) goes to zero and K1(x) goes to inf, and as x goes up, I1(x) goes to inf and K1(x) goes to zero. given that y(x) is finite, c1 and c2 go to inf or zero depending on their factor. and here the problem begins.

John D'Errico 2020-8-17

在 MATLAB Online 中打开

It is not that your followup question was boring. I am not able to watch answers every minute of the day however. I'll often only look in once a day, sometimes for a while, but then I will not return for some hours. We do have real lives to live, as you must recognize.

However, David is correct, in that your question is difficult to answer. Anyway, you have asked it separately now I see.

You ask about the choice of double precision arithmetic and a tool like VPA. As I showed, my own HPF is at leasrt similar in speed to the use of VPA. There will possibly be some cases where one may be faster, some slower. At least they are of the same order of magnitude, of that I am happy.

But to go to a tool like VPA or HPF is a huge jump. And there is no simple middle ground. You cannot easily make these tools faster. And sadly, when you do try to use those tools, you will likely find your computer is no longer easily able to make efficient use of the multiple cores on your computer. You should see that VPA computations are not automatically multi-threaded.

For example, this was large enough that it took some seconds to perform on my computer.

digits 100
A = vpa(rand(200));
B = inv(A);

Normally, any computation that will take any serious time to compute has at least a good chance that MATLAB will be able to multi-thread part of it. But that computation will not go beyond a single core. Symbolic computations are not automatically multi-threaded, even when they are large scale, at least not at this time.

It is almost always best to use good methods of numerical analysis to avoid the need to ever move into the realm of these high precision tools. If speed is a factor, just throwing high precision at the entire problem is a bad idea. At most, you may be able to gain by using such tools for a small piece that is not time consuuming. But large scale high precision computing is not easily achieved. Again, look to better numerical analysis methods to avoid the problems in the first place.

In the case of your differential equation, that may be to carefully transform the problem. If there is no possible way to achieve that, then you may need to recognize that big problems simply take big time or require bigger computers.

hosein Javan 2020-8-17

编辑：hosein Javan 2020-8-17

Bjorn Gustavsson. you're right. I'll give a more thorough explanation. the PDE is an inhomogeneous laplace equation in 2d plane:

where R(x) is a periodic function at x with its main period equal to "L". and at y1<y<y2.

the solution is:

where:

"L" is the (minimum) length of periodicity of f(x,y) which occurs in x coordinate.

"n" is the harmonic number.

the "

" is a particular solution.

the "

" and "

" are integral coefficients and are unknown. these are also y-dependant but in each area they are constants. as following:

by applying boundary conditions we can find the above integral coefficients. (it is quite complicated to explain here but in one sentence it concerns the function and its derivative continuity at each y0 to y6).

now after this we have "n" systems of linear equations in the form:

the elements of A matrix are one of like these:

the goal of the program is to:

compute the integral coefficients matrix [C_n]
compute the solution using the general solution and the integral coefficients [C_n]

the problem appears in both steps when computing

and the rest (the second kind and the zero-order bessel functions). especially at high harmonic numbers of "n", and high values of "p".

-------------------------------

as for unit scaling, I have already tried that. normally both x and y coordinates are in term of meters. now if we change the unit of y, then the unit of L scales the same as well and the result does not change. therefore the input arguement of the bessel functions is dimensionless. unless it is in terms of radians and maybe we can replace by something? not sure.
there's another thing I'm thinking that I'm not sure. we wanted the accuracy of f(x,y) at first place and not its bessel functions. maybe the accuracy of bessel functions and their coefficients when they multiply, cancel each other and result into a good accuracy of "f".

-------------------------------

however rigth now I can not set a value lower than L=0.5 when trying to solve for n=99 harmonics.

I hope this time I've explained exactly. sorry again for too much writing.

hosein Javan 2020-8-17

hello again Mr. D'Errico. yes you are right, my apologies. thanks for valuable information. as I said that speed is a goal itself I might not even think about symbolic or other high-precision data classes, I would certainly download your toolbox and give it a try, but for now I must try my best to remain at numeric realm as you said. there is however one thing I'm uncertain about. that is we want the precision of the solution, it may not be needed to have great accuracy of bessel functions and their coefficients. consider this example:

(1e300)/1e299

ans =

10

now we add the error at this high value wich is itself too huge: (eps(1e300) =1.487e+284)

>> (1e300+eps(1e300))/1e299

ans =

10

the result is very accurate. or at this small value:

(1e-300+eps(1e-300))/1e-299

ans =

0.1

this is however the double-precision capability. maybe the numeric variation of my function must be considered as well. because as I explained for Mr. Gustavsson, the resluts become quite unexpected when exceeding a certain value. if you can read it, it would fully explain.

the question is: given that f(x,y) requires a maximum error of

, how much accuracy is at least reqired for its solution form: bessel functions and their corresponding coefficients

at each "n"?

thanks again for spending time and effort,

Walter Roberson 2020-8-19

I wonder if the techniques in this paper would help:

https://www.researchgate.net/publication/230815104_Numerical_Calculation_of_Bessel_Functions

Bjorn Gustavsson 2020-8-20

Well not much more than the observation that if

when x goes to zero then that function does not work well in that interval. Perhaps it is enough to set the corresponding coefficients to zero and get a system of equations that solves the remaining coefficients.

请先登录，再进行评论。

Answer 2

Pavel Holoborodko 2020-8-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/580290-how-to-make-symbolic-computations-fast#answer_486629

编辑：Pavel Holoborodko 2020-8-29

在 MATLAB Online 中打开

I am developer of Advanpix toolbox. You mentioned it in the comments and that is why I decided to respond.

There are several misconceptions regarding computation of special functions (especially difficult ones, like Bessel or Hypergeometric or alike).

(a) First one is that such functions are considered easy to compute. You can frequently see the suggestions like - just use series, Pade, Chebyshebv or some other approximations.

(b) Second one is that well established software (NAG, MATLAB, Python, Julia, etc.) provide high-quality implementation for such functions.

Both are false. We have been working in the area for quite some time and published some of our findings, e.g. see [1,2]. Even double-precision software fails miseably to deliver proper accuracy.

When it comes to extended precision the situation gets even more difficult. No fixed degree approximation can be used, accuracy must be tracked on the-fly during the computations, cancellation effects must be detected, overflow/underflow must be handled properly, etc. etc. Combine this with different kinds of arguments and orders (complex, pure imaginary), principal values and branch cuts - and this can easily become the subject of solid research for several years.

Luckily you don't have to do all this, as we already did it for our toolbox. There is a reason why profesionally developed software costs money - it simply saves you ton of time by providing ready-to-use high-quality functionality with guaranteed accuracy.

Another advantage of our toolbox over "free" software is performance. We are several orders of magnitude faster than VPA, HPF, MP or even Maple, Mathematica, Julia and Python [3]. Multi-theading are natively built-in and used in all operations.

So, if you consider the quality and performance, then it becomes clear that there is no such thing as "free" software. Seemingly "free" software costs you a lot of time due to longer execution times, pain in setting up, etc. But most importanly your research might fail just because some "free" library provided sub-optimal accuracy in computations. That is pretty high cost.

Note 1. We didn't run speed comparison with HPF. But results are clear because of the author claim: "As I showed, my own HPF is at leasrt similar in speed to the use of VPA. There will possibly be some cases where one may be faster, some slower. At least they are of the same order of magnitude, of that I am happy." The VPA is the slowest extended precision software out there, we made comprehensive comparison with it [2].

Note 2. Using your example, this is how Advanpix toolbox computes it with 100 digits of accuracy on my computer:

>> mp.Digits(100);
>> tic;besseli(1,mp(100));toc;
Elapsed time is 0.000923 seconds.

Links.

[1] https://www.advanpix.com/2016/05/12/accuracy-of-bessel-functions-in-matlab/

[2] https://www.advanpix.com/documentation/users-manual/#Symbolic_Math_Toolbox_VPA_vs_Multiprecision_Computing_Toolbox

[3] (Japanese) https://qiita.com/hanaata/items/9a267bbf9f0547c807e7

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

how to make symbolic computations fast?

4 个评论
显示 2更早的评论隐藏 2更早的评论

回答（2 个）

14 个评论
显示 12更早的评论隐藏 12更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

how to make symbolic computations fast?

4 个评论 显示 2更早的评论隐藏 2更早的评论

回答（2 个）

14 个评论 显示 12更早的评论隐藏 12更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

14 个评论
显示 12更早的评论隐藏 12更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论