sum of products skipping an index

Question

Ioannis Andreou 2020-2-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/503142-sum-of-products-skipping-an-index

编辑： John D'Errico 2020-2-1

采纳的回答： David Hill

Hello everyone, I have difficulties efficiently implementing the following equation:

For Inputs: a,b,c 1xN doubles

I want to compute d 1xN double

such that

What makes it problematic is the skipping of the index in the product. Note, that computing the whole product and then dividing by the unwanted term does not work, since that would likely result in 0/0. I do have access to the symbolic math toolbox if that is of any help.

Thanks in advance :)

Ioannis

9 个评论
显示 7更早的评论隐藏 7更早的评论

John D'Errico 2020-2-1

在 MATLAB Online 中打开

That will make what I was thinking to do simple, but...

Next, as I look more carefully at what you are doing, I see that it looks suspiciously like a Largrange interpolating polynomial. Then I think about what you are asking to do. I guess its too early in the morning today, as I just realized that unless b and c are VERY carefully chosen vectors, this computation is meaningless in floating point arithmetic.

For example, suppose the average value of b(j) - c(i) is on the order of 2 or larger? Hmm, with N=10000, what is 2^10000?

2^10000
ans =
   Inf
2^2000
ans =
   Inf

So even for n on the order of a little over 1000, we will get an overflow into inf. That means, unless b(j) - c(i) is almost always approximately 1, you will get garbage. That means you almost surely CANNOT compute what you want to compute in double precision arithmetic.

Then again, I return to the question, is this something mathematically insane, like a Lagrange interpolating polynomial computation, where the polynomial must pass through ~10000 points? If that is your goal, then we need to back up agin. Why do you think you can do that? Why are you NOT using an intepolating spline?

Ioannis Andreou 2020-2-1

编辑：Ioannis Andreou 2020-2-1

Ok so these are the bounds of the input-values:

Also most elements in a and b will be 0. About c there's not much known.

John D'Errico 2020-2-1

编辑：John D'Errico 2020-2-1

在 MATLAB Online 中打开

The problem is b and c. a is irrelevant at the moment.

Are you saying that b lives in the open interval (-2,3)? And that c is completely unknown, that it can be anything? Actually, all that matters is the set of numbers for b, as that alone make what you are doing impossible in floating point arithmetic, at least directly. For example, consider this simple pair of vectors, of length only 10.

b = rand(1,10)*5 - 2
b =
       2.8154      0.73403      0.60568     -0.84203      0.44449       1.1203       1.3957    -0.022424     -0.16282       2.9399
>> c = linspace(-2,3,10)
c =
           -2      -1.4444     -0.88889     -0.33333      0.22222      0.77778       1.3333       1.8889       2.4444            3
bc = b' - c
bc =
       4.8154       4.2599       3.7043       3.1488       2.5932       2.0377       1.4821      0.92655        0.371     -0.18456
        2.734       2.1785       1.6229       1.0674      0.51181    -0.043749      -0.5993      -1.1549      -1.7104       -2.266
       2.6057       2.0501       1.4946      0.93901      0.38346      -0.1721     -0.72765      -1.2832      -1.8388      -2.3943
        1.158      0.60242     0.046861     -0.50869      -1.0643      -1.6198      -2.1754      -2.7309      -3.2865       -3.842
       2.4445       1.8889       1.3334      0.77782      0.22227     -0.33329     -0.88884      -1.4444           -2      -2.5555
       3.1203       2.5647       2.0092       1.4536      0.89808      0.34252     -0.21303     -0.76859      -1.3241      -1.8797
       3.3957       2.8401       2.2846        1.729       1.1735       0.6179     0.062344     -0.49321      -1.0488      -1.6043
       1.9776        1.422      0.86646      0.31091     -0.24465      -0.8002      -1.3558      -1.9113      -2.4669      -3.0224
       1.8372       1.2816      0.72607      0.17052     -0.38504     -0.94059      -1.4962      -2.0517      -2.6073      -3.1628
       4.9399       4.3844       3.8288       3.2732       2.7177       2.1621       1.6066        1.051      0.49547     -0.06009

So b is randomly chosen in the interval [-2,3). I chose c to be linear over the same interval.

The array b'-c is a subtraction table, an array that contains as the (i,j) element b(j) - c(i). Next, consider what the next line does:

bc = bc - diag(diag(bc)) + eye(10);

It replaces the diagonal elements with 1. Next, you want to compute that product, then multiply it each element with a corresponding value of a. But what is the product?

See that if N is much larger, on the order of 1e3 or 1e4, you will almost always get an overflow for the intermediate result. That is true as long as b ranges over the domain you claim it to live in.

You essentilly cannot compute what you want to compute, in double precision arithmetic, even though the scheme I have shown would allow you to do the desired computation for moderate sizes of N. 10000 will start to stretch things, because MATLAB will need to compute this intermediate array of size N^2. An array of that size will use 10000^2*8 bytes of memory, so roughly 800 megabytes. Doable on most computeres these days, but if you start to get much larger than that, you will quickly blow away the memory you have.

Now, I suppose you could do something as above, but using logs. Thus if

bc = b' - c;
bc = bc - diag(diag(bc)) + eye(10);
S = mod(sum(sign(bc) < 0,2),2);
bclog = sum(log(abs(bc)),2);
bclog
bclog =
       3.2063
      -1.2572
      0.30808
       1.4415
       2.3488
       2.1465
       3.1321
       0.4049
     0.024712
       7.1962

So bclog is the sum of the natual logs of those elements. I've saved the sign of the product in S, since that will cause problems with the logs. But again, if N is large, then you cannot do this computation, UNLESS you keep it in log form. It will virtually always cause overflows, as long as b has that large of a domain. If you then try to revert the log of the result by using exp, it will either underfow or overflow.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Hill 2020-2-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/503142-sum-of-products-skipping-an-index#answer_413261

在 MATLAB Online 中打开

You could try this:

d=zeros(1,N);
for i=1:N
    for k=1:N
        d(i)=d(i)+a(k)*prod(b([1:k-1,k+1:end])-c(i));
    end
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

sum of products skipping an index

9 个评论
显示 7更早的评论隐藏 7更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

sum of products skipping an index

9 个评论 显示 7更早的评论隐藏 7更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

9 个评论
显示 7更早的评论隐藏 7更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论