Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Summation and its condition

Subject: Summation and its condition

From: akihiro mizutani

Date: 16 Feb, 2009 09:10:03

Message: 1 of 13

Hello.

I’ve just started to use matlab.
And I’m stacked in the problem on summation condition.

Now I’m writing a code to get the result of following summation

Sum[(Aj*Ak-Bj*Bk)(Am*An-Bm*Bn)]

Where An and Bn are arrays respectively and length(An)=length(Bn)=N.
And the summation should be carried out only for those combinations of indices j, k, m, n that satisfy following conditions;

J+k=m+n, j<k, m<n, j<m.

At this, the indices may take any integer value but from 1 to N only.

Does anyone know how to define such conditions and carry out the summation for all combinations which satisfy the conditions?
I would be grateful if someone gave me a tip for solving the problem.

Thanks
Akihro

Subject: Summation and its condition

From: John D'Errico

Date: 16 Feb, 2009 10:42:02

Message: 2 of 13

"akihiro mizutani" <aki.mizutani@gmail.com> wrote in message <gnbahb$4id$1@fred.mathworks.com>...
> Hello.
>
> I’ve just started to use matlab.
> And I’m stacked in the problem on summation condition.
>
> Now I’m writing a code to get the result of following summation
>
> Sum[(Aj*Ak-Bj*Bk)(Am*An-Bm*Bn)]
>
> Where An and Bn are arrays respectively and length(An)=length(Bn)=N.
> And the summation should be carried out only for those combinations of indices j, k, m, n that satisfy following conditions;
>
> J+k=m+n, j<k, m<n, j<m.
>
> At this, the indices may take any integer value but from 1 to N only.
>
> Does anyone know how to define such conditions and carry out the summation for all combinations which satisfy the conditions?
> I would be grateful if someone gave me a tip for solving the problem.

If N is at all large, this will be difficult. It
involves O(N^4) terms in the combinatorial
sum, at least if I did my arithmetic correctly.

So how large is N?

John

Subject: Summation and its condition

From: akihiro mizutani

Date: 16 Feb, 2009 11:58:01

Message: 3 of 13

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gnbftq$m42$1@fred.mathworks.com>...
> "akihiro mizutani" <aki.mizutani@gmail.com> wrote in message <gnbahb$4id$1@fred.mathworks.com>...
> > Hello.
> >
> > I’ve just started to use matlab.
> > And I’m stacked in the problem on summation condition.
> >
> > Now I’m writing a code to get the result of following summation
> >
> > Sum[(Aj*Ak-Bj*Bk)(Am*An-Bm*Bn)]
> >
> > Where An and Bn are arrays respectively and length(An)=length(Bn)=N.
> > And the summation should be carried out only for those combinations of indices j, k, m, n that satisfy following conditions;
> >
> > J+k=m+n, j<k, m<n, j<m.
> >
> > At this, the indices may take any integer value but from 1 to N only.
> >
> > Does anyone know how to define such conditions and carry out the summation for all combinations which satisfy the conditions?
> > I would be grateful if someone gave me a tip for solving the problem.
>
> If N is at all large, this will be difficult. It
> involves O(N^4) terms in the combinatorial
> sum, at least if I did my arithmetic correctly.
>
> So how large is N?
>
> John

Thank you very much for a kind reply.
The number of N is up to 2024.

Even if it is difficult to solve the problem with such a big number, N,
I would be grateful if you show me your idea!

Thanks a lot.

Akihiro

Subject: Summation and its condition

From: John D'Errico

Date: 16 Feb, 2009 14:24:01

Message: 4 of 13

"akihiro mizutani" <aki.mizutani@gmail.com> wrote in message <gnbahb$4id$1@fred.mathworks.com>...
> Hello.
>
> I’ve just started to use matlab.
> And I’m stacked in the problem on summation condition.
>
> Now I’m writing a code to get the result of following summation
>
> Sum[(Aj*Ak-Bj*Bk)(Am*An-Bm*Bn)]
>
> Where An and Bn are arrays respectively and length(An)=length(Bn)=N.
> And the summation should be carried out only for those combinations of indices j, k, m, n that satisfy following conditions;
>
> J+k=m+n, j<k, m<n, j<m.
>
> At this, the indices may take any integer value but from 1 to N only.
>
> Does anyone know how to define such conditions and carry out the summation for all combinations which satisfy the conditions?
> I would be grateful if someone gave me a tip for solving the problem.
>
> Thanks
> Akihro

Ok, if you must solve this, I'll show you a way.
It is surely not the only way, nor even the best
way, but it will work reasonably for reasonable
problems. With N as large as 2000, it might fall
short, but that is a big number.

I'll pick some random vectors A and B to test
this on.

  N = 250;
  A = rand(1,N);
  B = rand(1,N);

Now for the algorithm. We will look at the
sum of j+k, fixing that sum to reduce your
problem, and then loop over the various
possible values of that index sum. Thus

% L = j + k = m + n
finalsum = 0;
for L = 2:(2*N)
  % L can be no smaller than 2, and no
  % larger than 2*N.
  % The possible values of j and k (as well
  % as m and n) given L = j + k
  j = 1:N;
  k = L - j;
  ind = (j>0) & (j<=N) & (k>0) & (k<=N) & (j < k);
  j = j(ind);
  k = k(ind);
  
  % now choose the corresponding m and n
  ind = 1:length(j);
  [indjk,indmn] = meshgrid(ind);
  indjk = indjk(:);
  indmn = indmn(:);
  
  % toss those out where j >= m
  ind = j(indjk) >= j(indmn);
  indjk(ind) = [];
  indmn(ind) = [];
  
  finalsum = finalsum + sum(...
     (A(j(indjk)).*A(k(indjk)) - B(j(indjk)).*B(k(indjk))).* ...
     (A(j(indmn)).*A(k(indmn)) - B(j(indmn)).*B(k(indmn))));
end

finalsum
finalsum =
    -0.065462

I don't know if this is correct, since I'd need
to carefully check the result and I don't have the
time right now. I did try it for N = 250, and the
time was reasonable. Even better, it looks like
I was off in my original O(N^4) estimate. This
should take something like O(2*N^3) operations
which is doable for N = 2024.

HTH,
John

Subject: Summation and its condition

From: Bruno Luong

Date: 16 Feb, 2009 20:12:02

Message: 5 of 13

"akihiro mizutani" <aki.mizutani@gmail.com> wrote in message <gnbahb$4id$1@fred.mathworks.com>...
> Hello.
>
> I’ve just started to use matlab.
> And I’m stacked in the problem on summation condition.
>
> Now I’m writing a code to get the result of following summation
>
> Sum[(Aj*Ak-Bj*Bk)(Am*An-Bm*Bn)]
>
> Where An and Bn are arrays respectively and length(An)=length(Bn)=N.
> And the summation should be carried out only for those combinations of indices j, k, m, n that satisfy following conditions;
>
> J+k=m+n, j<k, m<n, j<m.
>
> At this, the indices may take any integer value but from 1 to N only.
>
> Does anyone know how to define such conditions and carry out the summation for all combinations which satisfy the conditions?
> I would be grateful if someone gave me a tip for solving the problem.
>
> Thanks
> Akihro

Try this:

clear

% Set N to:
% <100 to check correctness
% 2024 to see it can work on your large array
N=100;

A=rand(N,1);
B=rand(N,1);

% Build a list of couple indices (i,j): 1<=i<j<=N
%
Aprod=A(:)*A(:).';
Bprod=B(:)*B(:).';
%
s=0;
for jpk=4:2*N-3
    AmnBmn=0;
    for j=floor((jpk-3)/2):-1:max(jpk-N,1)
        k=jpk-j;
        m = j+1;
        n = k-1;
        AmnBmn = AmnBmn + Aprod(m,n)-Bprod(m,n);
        s = s + (Aprod(j,k)-Bprod(j,k))*AmnBmn;
    end
end
fprintf('sum that can operate on big array=%g\n', s);

if N<=100
   % Brute force method
    n=N;
    idx=(1:n);
    [L J M]=ndgrid(2:2*n,idx,idx); % L is k+j = m+n
    K=L-J;
    N=L-M;
    good=J<K & M<N & J<M & K<=n;
    T=(A(J(good)).*A(K(good))-B(J(good)).*B(K(good))) .* ...
        (A(M(good)).*A(N(good))-B(M(good)).*B(N(good)));
    s2 = sum(T(:));
    fprintf('sum by bruteforce=%g\n', s);
    N=n;
end

% Bruno

Subject: Summation and its condition

From: Bruno Luong

Date: 16 Feb, 2009 20:26:02

Message: 6 of 13

% Hi have made a small improvement:

clear

% Set N to:
% <=100 to check correctness
% 2024 to see it can work on your large array
% I can go up to 7000
N=100;

A=rand(N,1);
B=rand(N,1);

AmB=A(:)*A(:).' - B(:)*B(:).';
s=0;
for jpk=4:2*N-3
    AmnBmn=0;
    for j=floor((jpk-3)/2):-1:max(jpk-N,1)
        k=jpk-j;
        AmnBmn = AmnBmn + AmB(j+1,k-1);
        s = s + AmB(j,k)*AmnBmn;
    end
end
fprintf('sum that can operate on big array=%g\n', s);

if N<=100
    n=N;
    idx=(1:n);
    [L J M]=ndgrid(2:2*n,idx,idx);
    K=L-J;
    N=L-M;
    good=J<K & M<N & J<M & K<=n;
    T=(A(J(good)).*A(K(good))-B(J(good)).*B(K(good))) .* ...
        (A(M(good)).*A(N(good))-B(M(good)).*B(N(good)));
    s2 = sum(T(:));
    fprintf('sum by bruteforce=%g\n', s);
    N=n;
end

% Bruno

Subject: Summation and its condition

From: Bruno Luong

Date: 16 Feb, 2009 20:46:03

Message: 7 of 13

Sorry !!!

Should read
> fprintf('sum by bruteforce=%g\n', s2); % s2 replaces s

% Bruno

Subject: Summation and its condition

From: Roger Stafford

Date: 16 Feb, 2009 21:52:02

Message: 8 of 13

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gnbsu1$5p9$1@fred.mathworks.com>...
> ......
> N = 250;
> A = rand(1,N);
> B = rand(1,N);
>
> Now for the algorithm. We will look at the
> sum of j+k, fixing that sum to reduce your
> problem, and then loop over the various
> possible values of that index sum. Thus
>
> % L = j + k = m + n
> finalsum = 0;
> for L = 2:(2*N)
> % L can be no smaller than 2, and no
> % larger than 2*N.
> % The possible values of j and k (as well
> % as m and n) given L = j + k
> j = 1:N;
> k = L - j;
> ind = (j>0) & (j<=N) & (k>0) & (k<=N) & (j < k);
> j = j(ind);
> k = k(ind);
>
> % now choose the corresponding m and n
> ind = 1:length(j);
> [indjk,indmn] = meshgrid(ind);
> indjk = indjk(:);
> indmn = indmn(:);
>
> % toss those out where j >= m
> ind = j(indjk) >= j(indmn);
> indjk(ind) = [];
> indmn(ind) = [];
>
> finalsum = finalsum + sum(...
> (A(j(indjk)).*A(k(indjk)) - B(j(indjk)).*B(k(indjk))).* ...
> (A(j(indmn)).*A(k(indmn)) - B(j(indmn)).*B(k(indmn))));
> end
>
> finalsum
> finalsum =
> -0.065462
>
> I don't know if this is correct, since I'd need
> to carefully check the result and I don't have the
> time right now. I did try it for N = 250, and the
> time was reasonable. Even better, it looks like
> I was off in my original O(N^4) estimate. This
> should take something like O(2*N^3) operations
> which is doable for N = 2024.
>
> HTH,
> John

  John, you'll be glad to learn that if you replace your summands with ones to get the total count of the number of j,k,m,n combinations satisfying those conditions, then the total I get with your code is in agreement with the formula

 ((2*N-3)*(2*N^2-6*N+1)+3*(-1)^N)/48,

which I obtained empirically but is almost certainly correct. It means the odds of your calculation being correct are heavily in your favor.

Roger Stafford

Subject: Summation and its condition

From: Bruno Luong

Date: 17 Feb, 2009 07:19:05

Message: 9 of 13

% Vectorize the inner loop :
% It takes about 10s to compute the sum for N=10000!

clear

N=10000;

A=rand(N,1);
B=rand(N,1);

tic
s=0;
for jpk=4:2*N-3
    j=floor((jpk-3)/2):-1:max(jpk-N,1);
    k=jpk-j;
    AmnBmn = A(j+1).*A(k-1) - B(j+1).*B(k-1);
    AjkBjk = A(j).*A(k) - B(j).*B(k);
    s = s + sum(AjkBjk.*cumsum(AmnBmn));
end
toc
fprintf('sum that can operate on big array=%g\n', s);

% Bruno

Subject: Summation and its condition

From: akihiro mizutani

Date: 17 Feb, 2009 09:16:03

Message: 10 of 13

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gndod9$stt$1@fred.mathworks.com>...
> % Vectorize the inner loop :
> % It takes about 10s to compute the sum for N=10000!
>
> clear
>
> N=10000;
>
> A=rand(N,1);
> B=rand(N,1);
>
> tic
> s=0;
> for jpk=4:2*N-3
> j=floor((jpk-3)/2):-1:max(jpk-N,1);
> k=jpk-j;
> AmnBmn = A(j+1).*A(k-1) - B(j+1).*B(k-1);
> AjkBjk = A(j).*A(k) - B(j).*B(k);
> s = s + sum(AjkBjk.*cumsum(AmnBmn));
> end
> toc
> fprintf('sum that can operate on big array=%g\n', s);
>
> % Bruno

Dear All,

Thank you for your help and efforts to improve them! I'm glad to be able to see various examples you posted. Your examples seem to derive correct answers as far as I have checked by myself.
I had not thought the summation could be computed with such a short computation time... super!
And now, I need to solve the other 4 conditions,
1. j=k+2n, k~=n
2. j+k=2n, j<k
3. j=3k
4. j+k+n=m, j<k<n

I'll try to solve them referring to your help but also, of course, any advice is always welcome!

Thank you very much!

Akihiro


 

Subject: Summation and its condition

From: John D'Errico

Date: 17 Feb, 2009 10:47:02

Message: 11 of 13

"akihiro mizutani" <aki.mizutani@gmail.com> wrote in message <gndv8j$k9j$1@fred.mathworks.com>...

> Dear All,
>
> Thank you for your help and efforts to improve them! I'm glad to be able to see various examples you posted. Your examples seem to derive correct answers as far as I have checked by myself.
> I had not thought the summation could be computed with such a short computation time... super!
> And now, I need to solve the other 4 conditions,
> 1. j=k+2n, k~=n
> 2. j+k=2n, j<k
> 3. j=3k
> 4. j+k+n=m, j<k<n
>
> I'll try to solve them referring to your help but also, of course, any advice is always welcome!
>

You have what is essentially a 4 dimensional problem,
that would normally require 4 nested loops, one for
each variable, {j, k, n, m}.

The trick that I provided was to use the equality
constraint

  j+k = n+m

to essentially eliminate one of those loops. This
reduces the problem to three dimensions. As such,
I had you loop on the value of (j+k). Inside that
loop, I used the additional information that if
j+k = n+m are both given as some number L,
then we can find k if we know j. This reduces
the problem one step further, down to two loops.

So for example, in one of your cases (the third
one shown), you have the constraint that j = 3*k.

If we know the value of k here, then we can compute
j directly. So this permits us to reduce the 4 loop
problem down to one with only 3 loops, on each of
k, n, and m.

Likewise, in the second case, you have j+k = 2*n.
Here if we know the value of n, then (j+k) is given
as 2*n. Now we must loop over three variables,
m, n, and one of either j or k.

Again, what we did was show you how to artfully
reduce the 4 explicit nested loops down to a
problem with fewer explicit loops, carefully using
your constraints.

John

Subject: Summation and its condition

From: Roger Stafford

Date: 17 Feb, 2009 16:48:02

Message: 12 of 13

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gndod9$stt$1@fred.mathworks.com>...
> % Vectorize the inner loop :
> % It takes about 10s to compute the sum for N=10000!
>
> clear
>
> N=10000;
>
> A=rand(N,1);
> B=rand(N,1);
>
> tic
> s=0;
> for jpk=4:2*N-3
> j=floor((jpk-3)/2):-1:max(jpk-N,1);
> k=jpk-j;
> AmnBmn = A(j+1).*A(k-1) - B(j+1).*B(k-1);
> AjkBjk = A(j).*A(k) - B(j).*B(k);
> s = s + sum(AjkBjk.*cumsum(AmnBmn));
> end
> toc
> fprintf('sum that can operate on big array=%g\n', s);
>
> % Bruno

  Well done Bruno! Both you and John have done well by Akihiro's problem. This morning when I awoke I had an idea of how the computation could be speeded up, only to find that you had already thought of it and more.

Roger Stafford

Subject: Summation and its condition

From: Bruno Luong

Date: 17 Feb, 2009 17:10:04

Message: 13 of 13

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gnepo2$18e$1@fred.mathworks.com>...

>
> Well done Bruno! Both you and John have done well by Akihiro's problem. This morning when I awoke I had an idea of how the computation could be speeded up, only to find that you had already thought of it and more.
>

Hi Roger and John,

Thanks! Knowing your capabilities, I have no doubt that you can come up with a similar solution soon or later. The idea is regrouping the sum is quite straightforward. The details in the implementation are a little harder to find. I have spent no less than half an hour yesterday to get it right. Fun problem.

Best,

Bruno

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us