Convolution of Two Vectors

[art297]

I don't see what you want.  When, in the calculation and updates, should A be applied and how?

[vioravis]

Hi Art,

I have given an example below. Please let me know if you need any other information. Thank you.

DATA test;

input X Y A Z;

cards;

1 10 0.1 0

2 20 0.2 0

3 30 0.3 0

4 40 0.4 0

5 50 0.5 0

;

run;

Iteration 1:

SUM(X[1:1]*Y[5–2+2:5]) = SUM(X[1:1]*Y[5:5]) = SUM([1]*[50]) = 50

X[2] = Z[2] + A[2] --> X[2] = 50 + 0.2 = 50.2

test dataset after iteration 1:

X     Y  A  Z

1    10 0.1 0

50.2 20 0.2 50

3    30 0.3 0

4    40 0.4 0

5    50 0.5 0

This test will be used in the second iteration as follows

Iteration 2:

SUM(X[1:2]*Y[5–3+2:5]) = SUM(X[1:2]*Y[4:5]) = SUM([1,50.2]*[40,50]) = SUM(40+2510) = 2550

X[3] = Z[3] + A[3] --> X[3] = 2550 + 0.3 = 2550.3

test dataset after iteration 2:

X      Y   A  Z

1      10 0.1 0

50.2   20 0.2 50

2550.3 30 0.3 2550

4      40 0.4 0

5      50 0.5 0

[art297 · Accepted solution]

The other responders are more versed than I am regarding the suitability and efficiency of the math involved.  However, that said, one minor correction you could make to your code in order to attain the results you want is:

data test;

  input X Y A Z;

  cards;

1 10 0.1 0

2 20 0.2 0

3 30 0.3 0

4 40 0.4 0

5 50 0.5 0

;

run;

data _null_;

  if 0 then set test nobs=nobs;

  CALL SYMPUT('NUMREC',nobs); /*** put # of records into NUMREC macro var ***/

  stop;                       /*** stop, got number of records ***/

run;

/* now read values of x and y into arrays, and then

   reread test and do calculations*/

data want (drop=_:);

    array _xval(&numrec); /* create three arrays with same

                                   number of elements as there are

                                   records in test

                          */

    array _yval(&numrec);

    array _Aval(&numrec);

    _i=0;

    do until (eof1);      /* load the array with a, x and y values */

      set test end=eof1;

      _i+1;

      _xval(_i)=x;

      _yval(_i)=y;

      _aval(_i)=a;

    end;

    _rec=-1;

    do until (eof2);      /* read in each record separately */

      set test end=eof2;

      _rec+1;

      z=0;

      do _i=1 to _rec;

        z+_xval(_i)*_yval(&numrec-_rec+_i);

      end;

      if _rec gt 0 then _xval(_rec+1)=sum(_aval(_rec+1),z);

      output;

    end;

run;

[Tom]

You might improve performance by doing the multiplication only once.

proc transpose data=test out=x prefix=x ;

   var x;

run;

proc transpose data=test out=y prefix=y;

   var y;

run;

data _null_;

   call symputx('NOBS',nobs);

   if 0 then set test nobs=nobs;

run;

data convolute;

   if _n_=1 then set x y;

   array x x1-x&nobs ;

   array y y1-y&nobs;

   array xy (&nobs,&nobs) _temporary_;

   if _n_=1 then do ;

        do i=1 to &nobs; do j=1 to &nobs ;

          xy(i,j)=x(i)*y(i);

        end; end;

   end;

***** now write summation code using indexes into xy 2-dimensional array ;

...

run;

[vioravis]

Thanks, Tom. I will try this out.

[vioravis]

Hi Tom,

The XY array you have defined holds the following matrix. Could you please explain what summation are you referring to on this XY array? Thank you.

        [ ,1]   [,2]   [,3]   [,4]   [,5]

[1,]   10     10    10    10    10

[2,]   40     40    40    40    40

[3,]   90     90    90    90    90

[4,]  160  160  160  160  160

[5,]  250  250  250  250  250

[Rick_SAS]

As much as I hate to be a naysayer, any approach that use the definition of the convolution will be many times slower than using a fast Fourier transform. From an efficiency point of view, these algorithms are not the best way to solve the problem. If you are going to be doing convolutions often, you might want to invest in SAS/IML.

I'm not very familar with SAS/ETS, but is there any procedure there that can be used to convolve vectors?

[FriedEgg]

I beleive there is a Spectral Analysis Procedure in ETS that can calculate FFT.

[Ksharp]

I do not agree with Rick Wicklin's opinion.

Array can offer you the fastest method.

If you are familiar with C language, then it will be more fast.

DATA test;
input X Y A Z;
cards;
1 10 0.1 0
2 20 0.2 0
3 30 0.3 0
4 40 0.4 0
5 50 0.5 0
;
run; 
%let dsid=%sysfunc(open(test));
%let nobs=%sysfunc(attrn(&dsid,nobs));
%let dsid=%sysfunc(close(&dsid));
data want(keep=x y a z);
 set test end=last;
 array xx{&nobs} _temporary_ ;
 array yy{&nobs} _temporary_ ;
 array aa{&nobs} _temporary_ ;
 array zz{&nobs} _temporary_ ;
 xx{_n_}=x;yy{_n_}=y;aa{_n_}=a;zz{_n_}=z;
 if last then do;
  do i=1 to &nobs-1; k=0; sum=0;
    do j=i to 1 by -1;
      k+1;
      sum=sum+xx{j}*yy{&nobs-k+1};
    end;
    zz{i+1}=sum; xx{i+1}=sum+aa{i+1};
  end;
  do l=1 to &nobs;
     x=xx{l};y=yy{l};a=aa{l};z=zz{l};
     output;
  end;
 end;
run;

 

Ksharp

[FriedEgg]

Traditional method requires at least N**2 calculations (so, for a 8,000 observation time series that is at least 64 million multiplications).  FFT algorithms computational load is more proportional to N*log2(N) (about 104,000 - a huge improvement).  I would have to agree with Rick.  The issue here is more that in SAS FFT is not available in products the users is licensed.  SAS/ETS and SAS/IML both have the ability to calculate FFT.  To incorporate FFT into base SAS one could possibly create a wrapper to a pre-complied C function (www.fftw.org for example) using the modulen function, it does require at least enough knowledge in C to setup the routine with the proper attribute table.

[Ksharp]

FriedEgg

That is the problem of algorithm, not the problem of coding.

If you can offer a better algorithm, Maybe I could code it too.

[vioravis]

Thanks to everyone who helped me with the solution.  I am currently using Art's solution but will be exploring Ksharp's solution as well in the near future.