sum of products skipping an index

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

[Ioannis Andreou]

So I tried doing it with a mix of funcionhandles and symbolic summation

prodSkip = @(v,k) prod([v(1:k-1), v(k+1:end)]); % product skipping k
syms ind
d = @(i) symsum(a(ind)*prodwok(b-c(i),ind), ind, 1, N));

Where then [d(1),..d(N)] would be the wanted result. However this doesn't work, as it produces the error:

Error using sym/subsindex (line 855)

Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be

symbolic variables, and function body must be sym expression.

Does anyone know how to avoid this?

[John D'Errico]

Very often, a simple solution may be good for smalll values of N, but impossible for N large. That is, if N is on the order of a few hundred or a thousand or so, then I can envision one solution. But, if N may be 1e5 or 1e6, the simple solution would become impossible. So before I even try to offer an answer, what is N? How large can it be?

Whenever you have a problem with an unknown size, please tell us what you expect that size to be, as otherwise it can make answering your problem a waste of time.

[Ioannis Andreou]

Overall the values of the inputted arrays depend on N too and my goal is to sample results for different values of N. So how big I will N let be is yet to be determined. But for a start assuming N < 1e4 would be reasonable.

[John D'Errico]

1e4 is starting to get a bit on the large side, but still doable. What Release are you using?

[Ioannis Andreou]

2019b

[John D'Errico]

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]

I am not doing interpolation or anything related to Lagrange-polynomial even though it does indeed look similiar.

About the size of the result your are right, I didn't think of that. The values of the arrays are (mostly) between -1 and 1

The formula above would only give an interim result which will then be divided by a number of approximately same magnitude, resulting in values of 'normal' size. However I forgot, that one might simply obtain exactly 0 as a result if N is too large, which would mean I would afterwards divide 0/0 which is not defined. So I assume I'll have to manage with N < 1e3. I will also provide more information about the input-values in a minute

[Ioannis Andreou]

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]

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.

Answer 1

[David Hill]

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