Multivariate Polynomials

From: Michael Friendly
Date: Mon, 7 Feb 1994 16:02:36 GMT
Newsgroups: comp.lang.apl
Subject: Multivariate polynomials- a puzzle

In regression applications, it is common to generate the columns of a polynomial in X to degree d by the operation

   X {null}.* 0,{iota}d

and these columns can be represented by the vector of powers, 0 1 2 ... d.

I have been stumped trying to generate a representation of the multivariate generalization of this, and hope there are some who enjoy a good APL puzzle.

Let v=number of variates, X1, X2, ... Xv, and
d=maximum degree term in the desired polynomial.

For v=2, a degree d=2 polynomial has terms X1, X1*2, X2, X2*2, X1xX2, which can be represented by the matrix below, whose columns correspond to X1, X2, and whose rows correspond to terms, the first of which is the intercept.

        0 0
        1 0
        2 0
        0 1
        0 2
        1 1

For v=2, d=3, there are additional rows

        3 0
        0 3
        1 2
        2 1

meshed in, where in general, all unique combinations which sum to <=d are included.

The puzzle is to generate such a matrix for arbitrary v and d. There ought to be an elegant APL solution for this problem, but I've been unable to find one. I have tried a direct approach, based on

        ( {null}.{times} ) / Z Z ... Z                  (d  Z's)

where Z = 1,X1,X2, ... ,[1.5]Xv as well, but this seems less fruitful than the indirect representation I described above.

(Please try to give replies in APL; I'm struggling too much with J.)


From: (Jim Weigang)
Newsgroups: comp.lang.apl
Date: 10 Feb 1994 03:23:07 GMT
Subject: Multivariate polynomial puzzle

I believe the integer matrix of powers can be obtained using the two statements:


(In index origin 1.) Building the matrix of variables raised to the different powers is not something I would be inclined to do with a one-liner. Such a solution is likely to involve lots of unnecessary arithmetic, be incomprehensible, or both. Here's a program that generates the matrix in a way that avoids unnecessary multiplication:

      {del}M{<-}D MVPOWERS X;A;C;I;J;K;M;N;P;Q;R;V
[1]    @Produces a multivariate power matrix
[2]    @ Arguments:
[3]    @     X  - data matrix.  Has a row for each observation and a column
[4]    @          for each variable.  Number of columns is assumed to be >0.
[5]    @     D  - degree of polynomial (assumed to be >0)
[6]    @ Result:
[7]    @     M  - matrix of powers.  Has a row for each row of X and columns
[8]    @          for the various powers and combinations of the variables.
[9]    @          For D=2 and 2=1{drop}{rho}X, the columns of the result are:
[10]   @             1, X[;1], X[;1]*2, X[;2], X[;1]{times}X[;2], X[;2]*2
[11]   @          The order of the columns is given by the matrix:
[12]   @             (D{>=}+/T){slashbar}T{<-}{transpose}(V{rho}D+1){+
[13]   @          This matrix has a row for each column of M, and a column
[14]   @          for each variable.  The numbers give the powers of each
[15]   @          variable.
[16]   @
[17]   V{<-}({rho}X)[2]           @ number of variables
[18]   N{<-}({rho}X)[1]           @ number of observations
[19]   M{<-}1,{times}\{transpose}(D,N){rho}X[;V] @ powers of the last variable
[20]   P{<-}{neg}1+{iota}1{drop}{rho}M           @ degree for each column of M
[21]   I{<-}V
[22]  L1:{->}(0=I{<-}I-1)/L4      @ Loop for each additional variable, last to first
[23]   R{<-}C{<-}X[;I]            @   the variable
[24]   A{<-}(1 0{times}{rho}M){rho}M @ accumulators for new data
[25]   Q{<-}{iota}0
[26]   J{<-}1                     @   current power of R
[27]  L2:                         @   Loop for each power of the variable
[28]   K{<-}1{drop}(D{>=}P+J)/{iota}1{drop}{rho}M {+
                               +} @     cols that will be combined with this power
[29]   A{<-}A,R,M[;K]{times}{transpose}(({rho}K),N){rho}R {+
                               +} @     additional columns
[30]   Q{<-}Q,J,P[K]+J            @     degree for the new cols
[31]   {->}(D<J{<-}J+1)/L3        @     quit if done
[32]   R{<-}R{times}C             @     next higher power
[33]   {->}L2                     @   Endloop
[34]  L3:M{<-}M,A                 @   append the new powers
[35]   P{<-}P,Q
[36]   {->}L1                     @ Endloop
[37]  L4:

After writing this, I realized how to write the nested version:


(The statement above is written for APL*PLUS II with EvLevel=1. {mix} is {disclose} with axis in APL2, and {split} is {enclose} with axis.)

Although I have an aversion to complicated nested array solutions like this, I must admit that the nested version took less time to write than the looping version. However, I wonder which version would be faster for another programmer to understand. For the comparison to be valid, the nested version should be written on multiple lines and comments added.

Timing tests show that the nested version runs faster when X has few rows, but that the looping version is faster when X has more than about 20 rows (for D=2 and V=3). As the number of rows gets large (say 1000), the speed ratio approaches a limit that depends on D and V. (The ratio is caused by unnecessary multiplications that are discarded by compression in the nested solution.) Some examples:

          D     V      Ratio  (nested:looping)
         ---   ---    ------
          2     2       1.5
          2     3       2.7
          2     4       5.4
          2     8     145.8
          4     3       3.6
          8     3       4.4
          4     8     789.1

For practical use, you'll want to use the looping version.


Home Page