Как вычислить двойную сумму в Matlab, где верхняя граница второй суммы является нижней границей первой суммы?

У меня большие трудности с реализацией функции в Matlab, которая вызывает другие функции, которые я запрограммировал в разных файлах .m. Часть, в которой я застрял, - это часть, где берется сумма по различным значениям, введенным в другую функцию, где также есть сумма внутри другой функции. Нижняя граница первой суммы является верхней границей второй суммы.

Функции: двойная сумма

У меня есть функция Hh (n, x), работающая правильно для n, введенного как вектор, и x, введенного как скаляр. Из-за векторизованного ввода n сумму по Hi внутри функции In можно быстро вычислить, вызвав sum(Hh(0:n,x)).

Я хочу сделать то же самое для функции In, но поскольку n внутри In теперь находится в диапазоне от 0 до k-1, а во внешней функции k находится в диапазоне от 1 до n, я не знаю, как вычислить эту двойную сумму, где внутренняя sum имеет нижнюю границу внешней суммы в качестве верхней границы. Я хочу вычислить эту двойную сумму как можно эффективнее, так как позже я хочу сделать много симуляций с этими формулами. Теперь я оцениваю функцию в n раз, сохраняя каждое значение в векторе, а затем беру сумму, которая требует больших вычислительных ресурсов...

Мой код Matlab для функции In:

function in = In(n,c,alphaa,betaa, delta)
ie = 0:n;
in = -(exp(alphaa*c)/alphaa)...
     .*sum((betaa/alphaa).^(n-ie).*Hh(ie,betaa*c-delta))...
     -(betaa/alphaa).^(n+1)
end

Код Matlab для внешней функции, назовем ее теперь функцией f:

function f = F(n,a,mu,sigma,eta1,T)
for k = 1: n
    vector(k) =  In(k-1,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T)));
end
f = sum(vector);
end

Как я могу сделать ввод n внутри In векторизованным, чтобы мне не нужно было хранить все введенные значения n отдельно, а затем вычислять сумму, а вычислять сумму непосредственно для введенного вектора n.

Любая помощь приветствуется, так как я серьезно застрял в данный момент! Заранее спасибо!


matlab performance sum vectorization memory-efficient
person Peter Lawrence    schedule 20.06.2017    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Итак, чтобы конкретно ответить на ваш вопрос (ввод n в векторе), я бы сделал следующее:

% 1. создать вектор ivect, идущий от 0 до k-1 для каждого k, и вектор kvect, соответствующий итерации ivect:

ivect=[];
kvect=[];
for k=1:n
    ivect = [ivect 0:k-1];
    kvect = [kvect (k-1)*ones(1,k)];
end

% 2. создайте функцию, которая находится внутри двух сумм, в зависимости от ваших итераций i и k:

    function in = In2(k,ie,c,alphaa,betaa, delta)
in = -(exp(alphaa*c)/alphaa)*((betaa/alphaa).^(k-ie).*func(ie,betaa*c-delta))-(betaa/alphaa).^(k+1)./(k+1);
end

% 3. тогда суммируем:

f=sum(In2(kvect,ivect,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T))))

Сначала я думал, что этот метод уменьшит вычислительное время при вычислении суммы (конечно пункт 1. затратный, но можно сделать один раз и использовать для нескольких симуляций без необходимости восстановления этих векторов), но даже сравнивая просто строки:

tic
for k = 1: n
    vector(k) =  In(k-1,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T)));
end
f = sum(vector)
toc

с

tic
f=sum(In2(kvect,ivect,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T))))
toc

Я получаю 1,029899 секунды для первого метода и 1,684432 секунды для второго, с n = 5000. Но, к счастью, я нашел почему! Это связано с тем, что этот метод на самом деле порождает больше арифметических манипуляций, так как в функции In2 я заставляю свой код вычислять -(betaa/alphaa).^(k+1)./(k+1) для каждого k и i, а этот член зависит только от k. Итак, теперь давайте определим функции In3 и In4 так, чтобы:

function in = In3(k,ie,c,alphaa,betaa, delta)
in = -(exp(alphaa*c)/alphaa)*((betaa/alphaa).^(k-ie).*Hh(ie,betaa*c-delta));    end

function in = In4(k,alphaa,betaa)
    in=-(betaa/alphaa).^(k+1);
end

и позвони

tic
f=sum(In3(kvect,ivect,a-mu*T,-eta1,-1/(sigma*sqrt(T)),-(sigma*eta1*sqrt(T)))) +sum(In4([1:n],-eta1,-1/(sigma*sqrt(T))))
toc

И вуаля, я получаю тот же результат за 0,929530 секунд для N=5000. Итак, подводя итог, вы выигрываете немного времени с этим последним методом, но, чтобы быть полностью честным, я бы сказал, что у вас есть два недостатка:

  1. Вычисление времени ivect и kvect в начале - но я настаиваю, что эти векторы могут быть построены один раз и навсегда, сохранены и использованы для нескольких симуляций -

  2. Стоимость памяти: ivect и kvect имеют размер n(n+1)/2, что много для больших n

Надеюсь, это помогло вам, и в любом случае спасибо за интересный вопрос, который, вероятно, послужит мне в моем собственном исследовании!

person Toghe    schedule 21.06.2017