how to make symbolic computations fast?

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?

[Bjorn Gustavsson]

You will get replies if you are more detailed and specific in your question.

[hosein Javan]

sorry at first I got no answers for a few hours so I thought my question was too long or confusing so I tried to shorten it. I hope this edit is clear enough, thanks.

[Steven Lord]

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]

Steven Lord. please read the comments. thanks.

Answer 1

[John D'Errico]

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.

[hosein Javan]

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.

[hosein Javan]

I see my question was too long and boring I assume. sorry.

[Bjorn Gustavsson]

I'd beg to differ, I still think your new question is till to vague. What are the PDEs? What are they supposed to model? What are the geometry and boundary-conditions? (Do you really need to calculate at spatial argument 1e-400? If x is in units of light-years, then x = 1e-400 ly corresponds to an awfully small distance...)

[John D'Errico]

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.

[Walter Roberson]

When array or vector computations are to be done, elementwise, or some matrix cases, then Yes, in theory parallel symbolic processing could be done. In practice, it is not usually done. Memory management turns out to be tricky for symbolic expressions.

[hosein Javan]

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:

1. compute the integral coefficients matrix [C_n]
2. 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]

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]

"maybe the accuracy of bessel functions and their coefficients when they multiply, cancel each other and result into a good accuracy of "f"."

Not in my experience. I find that when bessel functions individually get into large numbers, that you need higher precision for the terms in order to end up with final results that do not overflow and are not completely wrong. This probably implies that there would be a composite function that could more directly calculate the value without overflow, but it might be difficult to calculate.

[hosein Javan]

hello Mr. Roberson. the accuracy of bessel function calculated by "double" for its first 16 digits are fine, like in in this comparison between "vpa" that I stated at my main post. as for the accuracy of their coefficients I no longer get any warning saying that "rcond" is small, becasue I have fixed this matrix by scaling. however I still do not obtain accurate result for smaller "p" values.

Mr. Roberson ,what do you think about scaling 2*pi? can I convert radian to another larger-scale unit so that the input of bessel lowers? last night I was thinking about that "2*pi" is something trigonometric. however I doubt that in fourier series definition, "radian" can be convertible; since it is the ratio of circular line and straght line. what do you think we can define sin(1)=0 instead of sin(2*pi)=0 and therefor in fourier series 2*pi/T be converted to 1/T?

[Walter Roberson]

You spoke about intermediate values on the order of 1e-400, and your description suggests that you expect other values to be more or less on the order of 1e400 to balance out the product to roughly 1.

In my experience, when you start getting into numbers that large or that small, that you need extended precision, that 16 or so digits is not enough to get reasonable results even if the range is (somehow) extended enough to cover the 1e400 / 1e-400 values.

I have no idea whether scaling by 2*pi will work. I have a vague memory from many years ago that there are commonly two sets of spherical harmonic functions, one scaled by 2*pi relative to the other, but that was a long time ago and my memory on that detail is likely wrong.

[hosein Javan]

thanks for replying Mr. Roberson. maybe I should start thinking about that as Mr.John D'Errico said, bigger equations require greater calculations. what would you propose? should I do that or talk to a mathematician?

[hosein Javan]

any ideas Mr. Bjorn Gustavsson?

[Walter Roberson]

I wonder if the techniques in this paper would help:

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

[Bjorn Gustavsson]

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]

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