Your comments
Dear Mohsen,
Later today we will release new version of toolbox with even faster interp1. Now interp1 hot-spots are implemented directly using C++ and thus it is ~100 times faster overall compared to original version.
However, I am sure that the issue is not related speed, but rather to the properties of the method you use in your script. Is convergence guaranteed for all N?
For some N it converges very quickly, for others - no convergence even after thousands of iterations.
N=79:
>> interp1_test_mp
iteration = 1 crit = 6.966185e-01 time = 0.0154 sec
iteration = 2 crit = 1.328338e-01 time = 0.0167 sec
iteration = 3 crit = 4.058389e-01 time = 0.0207 sec
iteration = 4 crit = 9.628735e-01 time = 0.0136 sec
iteration = 5 crit = 7.685783e-02 time = 0.0138 sec
iteration = 6 crit = 5.120528e-02 time = 0.0138 sec
iteration = 7 crit = 6.838371e-02 time = 0.0135 sec
iteration = 8 crit = 4.207984e-02 time = 0.0134 sec
iteration = 9 crit = 4.638199e-02 time = 0.0281 sec
iteration = 10 crit = 2.107455e-02 time = 0.0167 sec
iteration = 11 crit = 1.724844e-02 time = 0.0173 sec
iteration = 12 crit = 5.805996e-03 time = 0.0172 sec
iteration = 13 crit = 2.771057e-03 time = 0.0164 sec
iteration = 14 crit = 1.432664e-03 time = 0.0149 sec
iteration = 15 crit = 1.115638e-03 time = 0.0156 sec
iteration = 16 crit = 1.431466e-03 time = 0.0173 sec
iteration = 17 crit = 1.182536e-03 time = 0.0148 sec
iteration = 18 crit = 7.638381e-04 time = 0.0145 sec
Cp_crit =
0.0007638381016874159474483705763685007
N=80:>> interp1_test_mp iteration = 1 crit = 6.966185e-01 time = 0.0136 sec iteration = 2 crit = 1.707076e-01 time = 0.0126 sec iteration = 3 crit = 1.819746e-01 time = 0.0126 sec .... iteration = 3644 crit = 1.867263e-02 time = 0.0130 sec iteration = 3645 crit = 1.840743e-02 time = 0.0134 sec ....Updated version will be available in a few hours.
Dear Mohsen,
Please download & install new version of toolbox: http://goo.gl/pMXV3
Now it has much faster interp1.
Also I have updated your script with small corrections to be more mp-fiendly:
rng('default');
% Grid length and useful vectors
N = 50;
kp = linspace(mp('1e-6'),mp(5),N)';
n = mp(ones(1,2));
% parameters & given values
gam = mp(2); beta = mp('0.96');
R = mp('1.05'); w0 = mp('3.44');
prob = mp('[0.9500, 0.0500; 0.6750, 0.3250]');
w = sort( w0*( mp('.995') + randn(N,1)/mp('7.5')));
W = [w, mp(zeros(N,1))];
% initialize
Cp = R*kp*n; k = Cp;
% convergence criterion
Cp_crit = 1;
% Interpolate until convergence
ii = 0;
while Cp_crit > 1e-3
start = tic;
Cp0 = Cp;
EMUp = (Cp.^(-gam))*prob';
C = (beta*R*EMUp).^(-1/gam);
k = (kp*n + C - W)/R;
% Use the relation between C and k to ropose a new vector Cp corresponding to kp
for i=1:2
% Update the function by interpolation.
Cp(:,i) = interp1(k(:,i), C(:,i), kp,'linear','extrap');
end
Cp_crit = max(max(abs(Cp0-Cp)./(1+abs(Cp))));
ii = ii+ 1;
stop = toc(start);
fprintf('iteration = %d\ttime = %8.4f sec\n',ii,stop);
end
disp(Cp_crit)
Older version has following timings for N=50:>> interp1_test_mp
iteration = 1 time = 0.9576 sec
iteration = 2 time = 0.9761 sec
iteration = 3 time = 1.0043 sec
iteration = 4 time = 1.0012 sec
iteration = 5 time = 1.0255 sec
iteration = 6 time = 1.0131 sec
iteration = 7 time = 1.0001 sec
iteration = 8 time = 0.9991 sec
iteration = 9 time = 1.0098 sec
iteration = 10 time = 1.0682 sec
iteration = 11 time = 1.0019 sec
iteration = 12 time = 1.0133 sec
iteration = 13 time = 1.0193 sec
iteration = 14 time = 1.0049 sec
iteration = 15 time = 1.0054 sec
iteration = 16 time = 1.0192 sec
iteration = 17 time = 1.0178 sec
iteration = 18 time = 1.0121 sec
Cp_crit =
0.0003433297868733290059017411707557392
New one is ~25 times faster:>> interp1_test_mp
iteration = 1 time = 0.0436 sec
iteration = 2 time = 0.0449 sec
iteration = 3 time = 0.0454 sec
iteration = 4 time = 0.0429 sec
iteration = 5 time = 0.0429 sec
iteration = 6 time = 0.0435 sec
iteration = 7 time = 0.0421 sec
iteration = 8 time = 0.0424 sec
iteration = 9 time = 0.0427 sec
iteration = 10 time = 0.0431 sec
iteration = 11 time = 0.0431 sec
iteration = 12 time = 0.0419 sec
iteration = 13 time = 0.0420 sec
iteration = 14 time = 0.0416 sec
iteration = 15 time = 0.0421 sec
iteration = 16 time = 0.0420 sec
iteration = 17 time = 0.0422 sec
iteration = 18 time = 0.0423 sec
Cp_crit =
0.0003433297868733290059017411707557392
But, still it might take a lot of time for N=200 to converge. Probably you need to re-design your algorithm to converge in less iterations.If you are Ok - I will make this thread public, it might be useful for other users.
Dear Mohsen,
You are absolutely right - tests show that interp1 is unacceptably slow!
The slowest code is index lookup of interval where x belongs while interpolation.
Now it is simple linear search and of course should be improved.
Basically we have to implement histcounts, which is not trivial.
Will try to improve this in near days and let you know about the result.
Updated version has been released, download link: http://goo.gl/pMXV3
Now I am working on speed optimization of Bessel functions for special cases.
For example new version already has much faster bessely(n,x) when n - integer, x - non-negative real in quadruple precision.
Please try it (if you use it in your work). I am finishing besselj - will be available very soon as well.
Bug has been fixed, updated version will be released in a few hours (after compilation finishes).
All relational operations were affected, and also element-wise power x.^y.
Dear Michael,
Thank you for your report, will fix it asap.
(Recently I have added support for sparse matrices in relational operations - probably something went wrong).
Just want to let you know, that I have extended toolbox with support of "shallow" complex numbers:
>> x = mp(complex(12))
x =
12 + 0i
>> isreal(x)
ans =
0But, in comparison to MATLAB, toolbox supports complex numbers much more correctly, including branch cuts, signed zeros as imaginary part, etc.
If you are interested in details, I wrote short post with examples for demonstration:
http://www.advanpix.com/2016/04/28/branch-cuts-and-signed-zeros-in-matlab/
Toolbox allows working with floating-point numbers from 9.5303e-323228497 to 5.2464e+323228495 with arbitrary precision (arbitrary length of mantissa):
>> format longG
>> mp.Digits(34);
>> mp('pi')
ans =
3.141592653589793238462643383279503
>> mp.Digits(100);
>> mp('pi')
ans =
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862
8034825342117068
>> mp.Digits(1000);
>> mp('pi')
ans =
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862
80348253421170679821480865132823066470938446095505822317253594081284811174502841027019385
21105559644622948954930381964428810975665933446128475648233786783165271201909145648566923
46034861045432664821339360726024914127372458700660631558817488152092096282925409171536436
78925903600113305305488204665213841469519415116094330572703657595919530921861173819326117
93105118548074462379962749567351885752724891227938183011949129833673362440656643086021394
94639522473719070217986094370277053921717629317675238467481846766940513200056812714526356
08277857713427577896091736371787214684409012249534301465495853710507922796892589235420199
56112129021960864034418159813629774771309960518707211349999998372978049951059731732816096
31859502445945534690830264252230825334468503526193118817101000313783875288658753320838142
06171776691473035982534904287554687311595628638823537875937519577818577805321712268066130
01927876611195909216420199
The range 9.5303e-323228497 to 5.2464e+323228495 is fixed but mantissa(precision) can be of any length.
The range covers pretty much all atoms in universe.
Could you show example where you need wider exponent range?
Yes, max/min numbers are the same for all precisions except IEEE quadruple.
Range of numbers depends on floating point exponent.
It is small for IEEE quadruple. For all other precisions we use the same large exponent.
Customer support service by UserEcho
Dear Hector,
Thank you very much for your report.
The error means that divide&conquer algorithm didn't converge.
Could you please share the input matrix so that I can reproduce the situation?
Thank you,
Pavel.