# sum of products skipping an index

16 views (last 30 days)
Ioannis Andreou on 1 Feb 2020
Edited: John D'Errico on 1 Feb 2020
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.
Ioannis
John D'Errico on 1 Feb 2020
Edited: John D'Errico on 1 Feb 2020
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?
prod(bc,2)
ans =
-24.687
-0.28446
-1.3608
4.2271
-10.473
8.5545
-22.922
-1.4991
-1.025
1334.4
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.

David Hill on 1 Feb 2020
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

### Categories

Find more on Number Theory in Help Center and File Exchange

R2019b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!