Dear Didier,

Thank you very much for the suggestion!

What would be the most important feature set?

(NFFT has quite a bit of features)

Integer matrices have to be converted to mp only if they participate in floating-point operations later on.

In this particular cases, just convert the matrices to mp after forming them as integers:

% 1. Compute elements of matrices as integers:
% ... (your code from above)

% 2. Convert them to mp:
A = mp(A);
G = mp(G);

% 3. Floating-point code with A and G goes afterwards:
% ...

P.S.

There is no integer matrices in MATLAB (unless you create them with explicit type - zeros(...,'int32')).

By default, zeros, ones, and everything else is created as double matrices. The floating-point is just not shown if number has no fractional part.

Dear Hector,

The A matrix contains NaN elements - and thus SVD cannot be computed:

>> mp.Digits(2000);
>> load Variables1.mat
>> isnan(A)
ans =
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     1     1     0     0
     0     0     1     0     1     1     1
     0     0     1     1     0     1     1
     0     0     0     1     1     0     0
     0     0     0     1     1     0     0

Next version of toolbox will show error message if there are NaN and Inf elements in a matrix.

Dear Hector,

Thank you for your assistance and report.

This issue has been fixed. Please use new version of toolbox - 4.0.0.11247.

Thank you.

Could you please send variables.mat to pavel@advanpix.com?

So that I will be able to reproduce the situation.

Yes, absolutely, this is very essential function. It can be called in two ways:

>> mp.Digits(34);
>> mp('eps')
ans = 
    1.925929944387235853055977942584927e-34
>> mp.eps
ans = 
    1.925929944387235853055977942584927e-34

There are also similar functions (and many others):

>> mp.realmax
ans = 
    1.189731495357231765085759326628007e+4932
>> mp.realmin
ans = 
    3.362103143112093506262677817321753e-4932

Please take a look to the last section on the page: http://www.advanpix.com/documentation/function-reference/

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.

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

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

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)

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

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

If you are Ok - I will make this thread public, it might be useful for other users.