4
\$\begingroup\$

I was trying to rewrite the Python code in MATLAB.

The result is consistent.

enter image description here

But, the MATLAB code is so slow. Any help would be appreciated.

Ref python code link

The MATLAB code written by me is as follows: so slow, whereas definitely correct.

I guess the critical time-costing step is compMultiFoxHIntegrand function. How can I rewrite it to make it faster?

function result = compMultiFoxH(params, nsubdivisions, boundaryTol)
    if nargin < 3
        boundaryTol = 0.0001;
    end

    boundaries = detBoundaries(params, boundaryTol);
    dim = numel(boundaries);
    
    signs = dec2bin(0:2^dim-1) - '0';
    signs(signs == 0) = -1;
    signs=signs.*(-1);

    inputs = repmat({(0:floor(nsubdivisions/2)-1)}, 1, dim);
    code = cartesian(inputs{:});
    
    quad = 0;
    % res = [];
    for i = 1:size(signs, 1)
        fprintf("i=%d\n",i);
        points = signs(i, :) .* (bsxfun(@plus, code, 0.5)) .* (boundaries * 2 / nsubdivisions);
        integrandVals = (compMultiFoxHIntegrand(points, params));
        % res = [res; integrandVals];
        quad = quad + sum(integrandVals);
    end
    volume = prod(2 * boundaries / nsubdivisions);
    result = quad * volume;
end

function boundaries = detBoundaries(params, tol)
    boundary_range = 0:0.05:50;
    dims = numel(params{1});
    boundaries = zeros(1, dims);
    for dim_l = 1:dims
        points = zeros(length(boundary_range), dims);
        points(:, dim_l) = boundary_range;
        abs_integrand = abs(compMultiFoxHIntegrand(points, params));
        index = find(abs_integrand > tol * abs_integrand(1), 1, 'last');
        boundaries(dim_l) = boundary_range(index);
    end
end

function result = compMultiFoxHIntegrand(y, params)
    z = params{1};
    mn = params{2};
    pq = params{3};
    c = params{4};
    d = params{5};
    a = params{6};
    b = params{7};

    m = cellfun(@(x) x(1), mn);
    n = cellfun(@(x) x(2), mn);
    p = cellfun(@(x) x(1), pq);
    q = cellfun(@(x) x(2), pq);
    npoints = size(y, 1);
    dims = size(y, 2);
    s = 1j * y;

    lower = zeros(1, dims);
    upper = zeros(1, dims);
    for dim_l = 1:dims
        if ~isempty(b{dim_l})
            bj = cellfun(@(x) x(1), b{dim_l}(1:m(dim_l+1)));
            Bj = cellfun(@(x) x(2), b{dim_l}(1:m(dim_l+1)));
            lower(dim_l) = -min(bj ./ Bj);
        else
            lower(dim_l) = -100;
        end
        if ~isempty(a{dim_l})
            aj = cellfun(@(x) x(1), a{dim_l}(1:n(dim_l+1)));
            Aj = cellfun(@(x) x(2), a{dim_l}(1:n(dim_l+1)));
            upper(dim_l) = min((1 - aj) ./ Aj);
        else
            upper(dim_l) = 1;
        end
    end
    mindist = norm(upper - lower);
    sigs = 0.5 * (upper + lower);
    for j = 1:n(1)
        num = 1 - c{j}(1) - sum(c{j}(2:end) .* lower);
        cnorm = norm(c{j}(2:end));
        newdist = abs(num) / (cnorm + eps);
        if newdist < mindist
            mindist = newdist;
            sigs = lower + 0.5 * num * c{j}(2:end) / (cnorm * cnorm);
        end
    end
    s = bsxfun(@plus, s, sigs);

    s1 = [ones(npoints, 1), s];
    prod_gam_num = 1 + 0j;
    prod_gam_denom = 1 + 0j;
    for j = 1:n(1)
        prod_gam_num = prod_gam_num .* gamma(sym(1 - s1 * c{j}'));
    end
    for j = 1:q(1)
        prod_gam_denom = prod_gam_denom .* gamma(sym(1 - s1 * d{j}'));
    end
    for j = n(1)+1:p(1)
        prod_gam_denom = prod_gam_denom .* gamma(sym(s1 * c{j}'));
    end
    for dim_l = 1:dims
        for j = 1:n(dim_l+1)
            prod_gam_num = prod_gam_num .* gamma(sym(1 - a{dim_l}{j}(1) - a{dim_l}{j}(2) * s(:, dim_l)));
        end
        for j = 1:m(dim_l+1)
            prod_gam_num = prod_gam_num .* gamma(sym(b{dim_l}{j}(1) + b{dim_l}{j}(2) * s(:, dim_l)));
        end
        for j = n(dim_l+1)+1:p(dim_l+1)
            prod_gam_denom = prod_gam_denom .* gamma(sym(a{dim_l}{j}(1) + a{dim_l}{j}(2) * s(:, dim_l)));
        end
        for j = m(dim_l+1)+1:q(dim_l+1)
            prod_gam_denom = prod_gam_denom .* gamma(sym(1 - b{dim_l}{j}(1) - b{dim_l}{j}(2) * s(:, dim_l)));
        end
    end

    zs = z .^ -s;
    result = double((prod_gam_num ./ prod_gam_denom) .* prod(zs, 2) ./ (2 * pi)^dims);
end



function C = cartesian(varargin)
    args = varargin;
    n = nargin;

    [F{1:n}] = ndgrid(args{:});

    for i=n:-1:1
        G(:,i) = F{i}(:);
    end

    C = unique(G , 'rows');
end

clear;clc;

% Example usage
params1 = {...
    [16.2982237081499, 16.2982237081499, 16.2982237081499, 16.2982237081499], ...
    {[0, 0], [2, 1], [2, 1], [2, 1], [2, 1]}, ...
    {[0, 1], [1, 2], [1, 2], [1, 2], [1, 2]}, ...
    {}, ...
    {[0, 1, 1, 1, 1]}, ...
    {{[1, 2]}, {[1, 2]}, {[1, 2]}, {[1, 2]}}, ...
    {{[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}} ...
};

result = compMultiFoxH(params1, 20);
format longG
disp(result);

enter image description here

\$\endgroup\$
0

2 Answers 2

Reset to default
5
\$\begingroup\$

There are some obvious inefficiencies in the code that I can point out.

The profiler shows that most time is spent doing symbolic math. You use symbolic math to compute the product of a series of gamma functions. There is no reason to compute this symbolically. MATLAB has a numeric gamma function. Simply leaving out the call to sym in your code should reduce the computational cost of your code by a few orders of magnitude.

Symbolic math is useful in a few circumstances. Mostly when doing symbolic manipulation, for example when finding the zeros of a function, which yields a new equation that can be evaluated many times over. The complicated part is done once, then the numerical evaluation that needs to be done millions of times becomes cheaper. In your case, all you do with the symbolic math is evaluate a numerical result, there is no symbolic manipulation. Numerical evaluation should be done numerically, it is a lot cheaper that way.

The other obvious inefficiency is all the unnecessary calls to cellfun. If you define the params data structure to be more convenient, then you won't need to copy and move around all these values every time you call the inner function compMultiFoxHIntegrand.

For example,

    mn = params{2};
    % ...
    m = cellfun(@(x) x(1), mn);
    n = cellfun(@(x) x(2), mn);

can be simplified by defining params{2} to be a 2D array instead of a cell array of 1D arrays. Your current implementation is an obvious 1-to-1 translation of the Python list of lists, but this is not very efficient at all. It is also not efficient in Python, the original code would also be more efficient using a 2D array.

    mn = params{2};
    m = mn(:, 1);
    n = mn(:, 2);
% ...

params1 = {...
    % ...
    [0, 0; 2, 1; 2, 1; 2, 1; 2, 1], ...
    % ...
};

Also, since you don't use mn in any other way but to separate it into its columns (at least not in the code you show), you could define params to have the two arrays separately. Note that m = mn(:, 1) means the data is being copied over.

If you use a struct instead of a cell array, this would be more convenient too:

params1 = [];
params1.m = [0, 2, 2, 2, 2];
params1.n = [0, 1, 1, 1, 1];

When specific indices into a cell array have specific meanings, it becomes very easy to mix things up. Making a change to the definition of params (e.g. by adding a value) could mean having to update all the code that uses it, and making mistakes there is very easy. A struct in MATLAB is very much like a cell array, but the indexing is done with a name instead of a number. This makes it much less likely that indices get mixed up, makes the code easier to read (e.g. you don't need to do mn = params{2}; just to give the value a readable name in the remainder of the code, you can just use params.mn everywhere, this indexing is very cheap), and makes it trivial to make changes to the definition of params.

\$\endgroup\$
2
\$\begingroup\$

Avoid Symbolic Math for Numerical Evaluations:

  • Replace symbolic computations with MATLAB's numeric gamma function. Symbolic math should only be used for symbolic manipulations, not for numerical evaluations.

It does work. Now it's much faster.

clear;clc;

% Example usage
params1 = {...
    [16.2982237081499, 16.2982237081499, 16.2982237081499, 16.2982237081499], ...
    {[0, 0], [2, 1], [2, 1], [2, 1], [2, 1]}, ...
    {[0, 1], [1, 2], [1, 2], [1, 2], [1, 2]}, ...
    {}, ...
    {[0, 1, 1, 1, 1]}, ...
    {{[1, 2]}, {[1, 2]}, {[1, 2]}, {[1, 2]}}, ...
    {{[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}} ...
};

result = compMultiFoxH(params1, 20);
format longG
disp(result);


function result = compMultiFoxH(params, nsubdivisions, boundaryTol)
    if nargin < 3
        boundaryTol = 0.0001;
    end

    boundaries = detBoundaries(params, boundaryTol);
    dim = numel(boundaries);
    
    signs = dec2bin(0:2^dim-1) - '0';
    signs(signs == 0) = -1;
    signs=signs.*(-1);

    inputs = repmat({(0:floor(nsubdivisions/2)-1)}, 1, dim);
    code = cartesian(inputs{:});
    
    quad = 0;
    % res = [];
    for i = 1:size(signs, 1)
        % fprintf("i=%d\n",i);
        points = signs(i, :) .* (bsxfun(@plus, code, 0.5)) .* (boundaries * 2 / nsubdivisions);
        integrandVals = (compMultiFoxHIntegrand(points, params));
        % res = [res; integrandVals];
        quad = quad + sum(integrandVals);
    end
    volume = prod(2 * boundaries / nsubdivisions);
    result = quad * volume;
end

function boundaries = detBoundaries(params, tol)
    boundary_range = 0:0.05:50;
    dims = numel(params{1});
    boundaries = zeros(1, dims);
    for dim_l = 1:dims
        points = zeros(length(boundary_range), dims);
        points(:, dim_l) = boundary_range;
        abs_integrand = abs(compMultiFoxHIntegrand(points, params));
        index = find(abs_integrand > tol * abs_integrand(1), 1, 'last');
        boundaries(dim_l) = boundary_range(index);
    end
end

function result = compMultiFoxHIntegrand(y, params)
    z = params{1};
    mn = params{2};
    pq = params{3};
    c = params{4};
    d = params{5};
    a = params{6};
    b = params{7};

    m = cellfun(@(x) x(1), mn);
    n = cellfun(@(x) x(2), mn);
    p = cellfun(@(x) x(1), pq);
    q = cellfun(@(x) x(2), pq);
    npoints = size(y, 1);
    dims = size(y, 2);
    s = 1j * y;

    lower = zeros(1, dims);
    upper = zeros(1, dims);
    for dim_l = 1:dims
        if ~isempty(b{dim_l})
            bj = cellfun(@(x) x(1), b{dim_l}(1:m(dim_l+1)));
            Bj = cellfun(@(x) x(2), b{dim_l}(1:m(dim_l+1)));
            lower(dim_l) = -min(bj ./ Bj);
        else
            lower(dim_l) = -100;
        end
        if ~isempty(a{dim_l})
            aj = cellfun(@(x) x(1), a{dim_l}(1:n(dim_l+1)));
            Aj = cellfun(@(x) x(2), a{dim_l}(1:n(dim_l+1)));
            upper(dim_l) = min((1 - aj) ./ Aj);
        else
            upper(dim_l) = 1;
        end
    end
    mindist = norm(upper - lower);
    sigs = 0.5 * (upper + lower);
    for j = 1:n(1)
        num = 1 - c{j}(1) - sum(c{j}(2:end) .* lower);
        cnorm = norm(c{j}(2:end));
        newdist = abs(num) / (cnorm + eps);
        if newdist < mindist
            mindist = newdist;
            sigs = lower + 0.5 * num * c{j}(2:end) / (cnorm * cnorm);
        end
    end
    s = bsxfun(@plus, s, sigs);

    s1 = [ones(npoints, 1), s];
    prod_gam_num = 1 + 0j;
    prod_gam_denom = 1 + 0j;
    for j = 1:n(1)
        prod_gam_num = prod_gam_num .* gammac( (1 - s1 * c{j}'));
    end
    for j = 1:q(1)
        prod_gam_denom = prod_gam_denom .* gammac( (1 - s1 * d{j}'));
    end
    for j = n(1)+1:p(1)
        prod_gam_denom = prod_gam_denom .* gammac( (s1 * c{j}'));
    end
    for dim_l = 1:dims
        for j = 1:n(dim_l+1)
            prod_gam_num = prod_gam_num .* gammac( (1 - a{dim_l}{j}(1) - a{dim_l}{j}(2) * s(:, dim_l)));
        end
        for j = 1:m(dim_l+1)
            prod_gam_num = prod_gam_num .* gammac( (b{dim_l}{j}(1) + b{dim_l}{j}(2) * s(:, dim_l)));
        end
        for j = n(dim_l+1)+1:p(dim_l+1)
            prod_gam_denom = prod_gam_denom .* gammac( (a{dim_l}{j}(1) + a{dim_l}{j}(2) * s(:, dim_l)));
        end
        for j = m(dim_l+1)+1:q(dim_l+1)
            prod_gam_denom = prod_gam_denom .* gammac( (1 - b{dim_l}{j}(1) - b{dim_l}{j}(2) * s(:, dim_l)));
        end
    end

    zs = z .^ -s;
    result = double((prod_gam_num ./ prod_gam_denom) .* prod(zs, 2) ./ (2 * pi)^dims);
end



function C = cartesian(varargin)
    args = varargin;
    n = nargin;

    [F{1:n}] = ndgrid(args{:});

    for i=n:-1:1
        G(:,i) = F{i}(:);
    end

    C = unique(G , 'rows');
end

function [f] = gammac(z)
% GAMMA  Gamma function valid in the entire complex plane.
%        Accuracy is 15 significant digits along the real axis
%        and 13 significant digits elsewhere.
%        This routine uses a superb Lanczos series
%        approximation for the complex Gamma function.
%
%        z may be complex and of any size.
%        Also  n! = prod(1:n) = gamma(n+1)
%
%usage: [f] = gamma(z)
%
%tested on versions 6.0 and 5.3.1 under Sun Solaris 5.5.1
%
%References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
%            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
%            Y. Luke, "Algorithms ... functions", 1977
%            J. Spouge,  SIAM JNA 31, 1994. pp. 931-944
%            W. Press,  "Numerical Recipes"
%            S. Chang, "Computation of special functions", 1996
%            W. J. Cody "An Overview of Software Development for Special
%            Functions", 1975
%
%see also:   GAMMA GAMMALN GAMMAINC PSI
%see also:   mhelp GAMMA
%
%Paul Godfrey
%[email protected]
%http://winnie.fit.edu/~gabdo/gamma.txt
%Sept 11, 2001
siz = size(z);
z=z(:);
zz=z;
f = 0.*z; % reserve space in advance
p=find(real(z)<0);
if ~isempty(p)
    z(p)=-z(p);
end
% 15 sig. digits for 0<=real(z)<=171
% coeffs should sum to about g*g/2+23/24
g=607/128; % best results when 4<=g<=5
c = [  0.99999999999999709182;
    57.156235665862923517;
    -59.597960355475491248;
    14.136097974741747174;
    -0.49191381609762019978;
    .33994649984811888699e-4;
    .46523628927048575665e-4;
    -.98374475304879564677e-4;
    .15808870322491248884e-3;
    -.21026444172410488319e-3;
    .21743961811521264320e-3;
    -.16431810653676389022e-3;
    .84418223983852743293e-4;
    -.26190838401581408670e-4;
    .36899182659531622704e-5];
%Num Recipes used g=5 with 7 terms
%for a less effective approximation
z=z-1;
zh =z+0.5;
zgh=zh+g;
%trick for avoiding FP overflow above z=141
zp=zgh.^(zh*0.5);
ss=0.0;
for pp=size(c,1)-1:-1:1
    ss=ss+c(pp+1)./(z+pp);
end
%sqrt(2Pi)
sq2pi=  2.5066282746310005024157652848110;
f=(sq2pi*(c(1)+ss)).*((zp.*exp(-zgh)).*zp);
f(z==0 | z==1) = 1.0;
%adjust for negative real parts
if ~isempty(p)
    f(p)=-pi./(zz(p).*f(p).*sin(pi*zz(p)));
end
%adjust for negative poles
p=find(round(zz)==zz & imag(zz)==0 & real(zz)<=0);
if ~isempty(p)
    f(p)=Inf;
end
f=reshape(f,siz);
return
endfunction [outputArg1,outputArg2] = untitled2(inputArg1,inputArg2)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
outputArg1 = inputArg1;
outputArg2 = inputArg2;
end
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.