Thread Subject: Vectorization for Quadratic Polynomial Regression

I have an original predictor matrix X0 with

size(X0) = [ m n0]

I'd like to perform the

linear-coefficient/quadratic-variable regression

y = X*W + e

size(X) = [ m n]

n = n0*(n0+3)/2

and X contains the columns of X0,
their squares, and their cross products.

Currently, X is created by the following double loop

X = X0;
for i = 1:n0-1
    for j= i+1:n0
        X = [ X X0(:,i).*X0(:,j)];
    end
end
X = [X X0.^2];

Example:

X0 = [ 1 2 3; 4 5 6]

X =
     1 2 3 2 3 6 1 4 9
     4 5 6 20 24 30 16 25 36

Can the double loop be vectorized?

TIA,

Greg

On Jul 1, 3:05 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> I have an original predictor matrix X0 with
>
> size(X0) = [ m n0]
>
> I'd like to perform the
>
> linear-coefficient/quadratic-variable regression
>
> y = X*W + e
>
> size(X) = [ m n]
>
> n = n0*(n0+3)/2
>
> and X contains the columns of X0,
> their squares, and their cross products.
>
> Currently, X is created by the following double loop
>
> X = X0;
> for i = 1:n0-1
> for j= i+1:n0
> X = [ X X0(:,i).*X0(:,j)];
> end
> end

X = [ ones(m,1) X X0.^2 ];

> Example:
>
> X0 = [ 1 2 3; 4 5 6]
>
> X =

  1 1 2 3 2 3 6 1 4 9
  1 4 5 6 20 24 30 16 25 36


> Can the double loop be vectorized?

TIA,

Greg

>> Can the double loop be vectorized?

Greg, I offer this not as an example of good programming, but of something
that works and that might give you some ideas:

terms = fullfact([3 3 3])-1;
terms(sum(terms,2)>2,:) = []
x=[ 1 2 3; 4 5 6]
exp(log(x)*terms')

-- Tom

On Jul 1, 4:32 pm, "Tom Lane" <tl...@mathworks.com> wrote:
> >> Can the double loop be vectorized?
>
> Greg, I offer this not as an example of good programming, but of something
> that works and that might give you some ideas:
>
> terms = fullfact([3 3 3])-1;
> terms(sum(terms,2)>2,:) = []
> x=[ 1 2 3; 4 5 6]
> exp(log(x)*terms')

This is ~ 4.5 times slower than the double loop
for 10,000 trials using the above example.

How is this generalized for

size(X0) = [ m n0] ?


X0 = repmat( [ 1 2 3 ; 4 5 6 ], 2, 2 );
[m n0] = size(X0) % [4 6]

tic
for k = 1:10000
    X = X0;
    for i = 1:n0-1
        for j= i+1:n0
            X = [ X X0(:,i).*X0(:,j)];
        end
    end
    X = [ones(m,1) X X0.^2];
end
[m n] = size(X) % [ 4 28 ]
toc % ~ 0.78 sec

tic
for k=1:1000
    terms = fullfact( [ ? ? ? ] ) - 1;
    terms( sum( terms, 2 ) > 2, : ) = [];
     Y = exp( log (X0 ) * terms' );
end
toc % ???

TIA,

Greg

> How is this generalized for
>
> size(X0) = [ m n0] ?

Greg, you can replace: [? ? ?]
by this: 3*ones(1,n0)

The idea is that we're getting all combinations of [0,1,2] as powers of the
x columns, and weeding out the rows where the sum of the exponents is larger
than the desired value of 2.

Again, I offer this as a curiosity rather than a real recommendation. For
one thing, it's not going to do good things if there are any negative values
in X0.

-- Tom

"Tom Lane" <tlane@mathworks.com> wrote in message
<g4ennn$l7t$1@fred.mathworks.com>...
> > How is this generalized for
> >
> > size(X0) = [ m n0] ?
>
> Greg, you can replace: [? ? ?]
> by this: 3*ones(1,n0)
>
> The idea is that we're getting all combinations of [0,1,2] as powers of the
> x columns, and weeding out the rows where the sum of the exponents is
larger
> than the desired value of 2.
>
> Again, I offer this as a curiosity rather than a real recommendation. For
> one thing, it's not going to do good things if there are any negative values
> in X0.

For gods sake, why do it the hard way?

Use bsxfun around realpow.

It won't care then even if you do have negative
numbers.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message
<g4ep5e$8f0$1@fred.mathworks.com>...
> "Tom Lane" <tlane@mathworks.com> wrote in message
> <g4ennn$l7t$1@fred.mathworks.com>...
> > > How is this generalized for
> > >
> > > size(X0) = [ m n0] ?
> >
> > Greg, you can replace: [? ? ?]
> > by this: 3*ones(1,n0)
> >
> > The idea is that we're getting all combinations of [0,1,2] as powers of the
> > x columns, and weeding out the rows where the sum of the exponents is
> larger
> > than the desired value of 2.
> >
> > Again, I offer this as a curiosity rather than a real recommendation. For
> > one thing, it's not going to do good things if there are any negative
values
> > in X0.
>
> For gods sake, why do it the hard way?
>
> Use bsxfun around realpow.
>
> It won't care then even if you do have negative
> numbers.
>
> John

Sorry. power is adequate here, and is
supported by bsxfun.

bsxfun(@power,(-1:.1:1)',0:5)

John

On Jul 1, 10:04=A0pm, "Tom Lane" <tl...@mathworks.com> wrote:
> > How is this generalized for
>
> > size(X0) =3D [ m n0] ?
>
> Greg, you can replace: =A0[? ? ?]
> by this: =A03*ones(1,n0)

Thanks. Now getting ~ 0.47 sec vs 0.78 sec for the loop version.

> The idea is that we're getting all combinations of [0,1,2] as powers of th=
e
> x columns, and weeding out the rows where the sum of the exponents is larg=
er
> than the desired value of 2.

OK. ... I'll have to come back to better understand fullfact.

> Again, I offer this as a curiosity rather than a real recommendation. =A0F=
or
> one thing, it's not going to do good things if there are any negative valu=
es
> in X0.

My original values are all positive. However, I will probably
standardize.

Greg

> For gods sake, why do it the hard way?
>
> Use bsxfun around realpow.

John, it's because I was using matrix multiplication rather than
element-by-element multiplication. This allowed me in effect to
exponentiate and multiply in one operation without creating a giant array.

Again I will say I offer this as a curiosity. There are many things wrong
with it, starting with these:

1. Negative x values will probably lead to a complex result with a tiny
imaginary part.

2. I created a potentially giant terms matrix only to weed out most of its
rows.

3. I am needlessly including variables to the 0th power in every term.

Just hoping it will lead Greg down a fruitful path, or that someone else
will improve upon it.

-- Tom

Greg Heath <heath@alumni.brown.edu> wrote in message
<cddc4bd3-1f78-4fad-90d2-a7cde7ec5189@e39g2000hsf.googlegroups.com>...
> I have an original predictor matrix X0 with
>
> size(X0) = [ m n0]
>
> I'd like to perform the
>
> linear-coefficient/quadratic-variable regression
>
> y = X*W + e
>
> size(X) = [ m n]
>
> n = n0*(n0+3)/2
>
> and X contains the columns of X0,
> their squares, and their cross products.
>
> Currently, X is created by the following double loop
>
> X = X0;
> for i = 1:n0-1
> for j= i+1:n0
> X = [ X X0(:,i).*X0(:,j)];
> end
> end
> X = [X X0.^2];
>
> Example:
>
> X0 = [ 1 2 3; 4 5 6]
>
> X =
> 1 2 3 2 3 6 1 4 9
> 4 5 6 20 24 30 16 25 36
>
> Can the double loop be vectorized?
>
> TIA,
>
> Greg
Hi,

It seems that you have messed up the formulation of the
quadratic least square problem. Isnt it(?):

yi=c0+c1*xi+c2*xi.^2

A=[1+0*x(:) x(:) x(:).^2];
c=y\A

%\= least square solution (check help in matlab)

On Jul 2, 2:16=A0pm, "Per Sundqvist" <sun...@fy.chalmers.se> wrote:
> Greg Heath <he...@alumni.brown.edu> wrote in message
>
> <cddc4bd3-1f78-4fad-90d2-a7cde7ec5...@e39g2000hsf.googlegroups.com>...
>
>
>
> > I have an original predictor matrix X0 with
>
> > size(X0) =3D [ m n0]
>
> > I'd like to perform the
>
> > linear-coefficient/quadratic-variable regression
>
> > y =3D X*W + e
>
> > size(X) =3D [ m n]
>
> > n =3D n0*(n0+3)/2
>
> > and X contains the columns of X0,
> > their squares, and their cross products.
>
> > Currently, X is created by the following double loop
>
> > X =3D X0;
> > for i =3D 1:n0-1
> > =A0 =A0 for j=3D i+1:n0
> > =A0 =A0 =A0 =A0 X =3D [ X X0(:,i).*X0(:,j)];
> > =A0 =A0 end
> > end
> > X =3D [X X0.^2];
>
> > Example:
>
> > X0 =3D [ =A01 2 3; 4 5 6]
>
> > X =3D
> > =A0 =A0 =A01 =A0 =A0 2 =A0 =A0 3 =A0 =A0 2 =A0 =A0 3 =A0 =A0 6 =A0 =A0 1=
 =A0 =A0 4 =A0 =A0 9
> > =A0 =A0 =A04 =A0 =A0 5 =A0 =A0 6 =A0 =A020 =A0 =A024 =A0 =A030 =A0 =A016=
 =A0 =A025 =A0 =A036
>
> > Can the double loop be vectorized?
>
> > TIA,
>
> > Greg
>
> Hi,
>
> It seems that you have messed up the formulation of the
> quadratic least square problem. Isnt it(?):
>
> yi=3Dc0+c1*xi+c2*xi.^2

No.

There are multiple variables ( e.g., n0 =3D 3 ) in X0, yielding
n0 original variable columns, n0 variable squared columns
and n0*(n0-1)/2 crossproduct columns in X.

When n0 =3D 3, n =3D 9 and, for size(y) =3D [m 1]
size(W) =3D [n+1 1] =3D [10 1]. For measurement i

yhati =3D W1*1 + W2*x1i + W3*x2i + W4*x3i +
              W5*x1i*x2i + W6 *x1i*x3i + W7*x2i*x3i +
              W8*x1i^2 + W9*x2i^2 + W10*x3i^2

 where

W =3D [ones(m,1) X]\y;

Hope this helps.

Greg

On Jul 2, 2:16=A0pm, "Per Sundqvist" <sun...@fy.chalmers.se> wrote:
> Greg Heath <he...@alumni.brown.edu> wrote in message
>
> <cddc4bd3-1f78-4fad-90d2-a7cde7ec5...@e39g2000hsf.googlegroups.com>...
>
>
>
> > I have an original predictor matrix X0 with
>
> > size(X0) =3D [ m n0]
>
> > I'd like to perform the
>
> > linear-coefficient/quadratic-variable regression
>
> > y =3D X*W + e
>
> > size(X) =3D [ m n]
>
> > n =3D n0*(n0+3)/2
>
> > and X contains the columns of X0,
> > their squares, and their cross products.
>
> > Currently, X is created by the following double loop
>
> > X =3D X0;
> > for i =3D 1:n0-1
> > =A0 =A0 for j=3D i+1:n0
> > =A0 =A0 =A0 =A0 X =3D [ X X0(:,i).*X0(:,j)];
> > =A0 =A0 end
> > end
> > X =3D [X X0.^2];
>
> > Example:
>
> > X0 =3D [ =A01 2 3; 4 5 6]
>
> > X =3D
> > =A0 =A0 =A01 =A0 =A0 2 =A0 =A0 3 =A0 =A0 2 =A0 =A0 3 =A0 =A0 6 =A0 =A0 1=
 =A0 =A0 4 =A0 =A0 9
> > =A0 =A0 =A04 =A0 =A0 5 =A0 =A0 6 =A0 =A020 =A0 =A024 =A0 =A030 =A0 =A016=
 =A0 =A025 =A0 =A036
>
> > Can the double loop be vectorized?
>
> > TIA,
>
> > Greg
>
> Hi,
>
> It seems that you have messed up the formulation of the
> quadratic least square problem.

Oh, I see you were misled by the use of the term
"Quadratic" instead of "Second Order" .

Learn from this:

Don't assume titles are precise indicators of thread content

(;>).

Greg

On Jul 2, 12:07=A0pm, "Tom Lane" <tl...@mathworks.com> wrote:
> > For gods sake, why do it the hard way?
>
> > Use bsxfun around realpow.
>
> John, it's because I was using matrix multiplication rather than
> element-by-element multiplication. =A0This allowed me in effect to
> exponentiate and multiply in one operation without creating a giant array=
.
>
> Again I will say I offer this as a curiosity. =A0There are many things wr=
ong
> with it, starting with these:
>
> 1. =A0Negative x values will probably lead to a complex result with a tin=
y
> imaginary part.
>
> 2. =A0I created a potentially giant terms matrix only to weed out most of=
 its
> rows.
>
> 3. =A0I am needlessly including variables to the 0th power in every term.
>
> Just hoping it will lead Greg down a fruitful path, or that someone else
> will improve upon it.

Thanks Guys.

 Although I can't use Tom's solution because of standardization,
I can't use John's either because bxsfun is not available in ver
6.5.1.

Since my problem only has n0 =3D 8 predictors, using the double loop
to get n =3D 44 is not affecting the model development progress.
However, not having a vectorized code is affecting my sleep!

Greg

>> Again I will say I offer this as a curiosity. There are many things wrong
>> with it, starting with these:
...
> However, not having a vectorized code is affecting my sleep!

Greg, I can't handle the stigma of my lousy solution and my concern over
your lack of sleep. Try this:

function terms=interactionterms(x)
[n,p] = size(x);
m = repmat(1:p,p,1); % form all combinations
mt = m'; % second half of each pair
lower = ~triu(ones(p)); % get only (i,j) with i<j
i = m(lower);
j = mt(lower);
terms = [ones(n,1), ... % constant term
         x, ... % linear terms
         x.^2, ... % squared terms
         x(:,i).*x(:,j)]; % pairwise interactions

This won't generalize to higher powers and products like my lousy solution,
but it's better for your problem.

-- Tom

Greg Heath <heath@alumni.brown.edu> wrote in message
<6df4a42e-3d83-4dd7-9441-c371361dab06@c65g2000hsa.googlegroups.com>...
> On Jul 2, 2:16=A0pm, "Per Sundqvist"
<sun...@fy.chalmers.se> wrote:
> > Greg Heath <he...@alumni.brown.edu> wrote in message
> >
> >
<cddc4bd3-1f78-4fad-90d2-a7cde7ec5...@e39g2000hsf.googlegroups.com>...
> >
> >
> >
> > > I have an original predictor matrix X0 with
> >
> > > size(X0) =3D [ m n0]
> >
> > > I'd like to perform the
> >
> > > linear-coefficient/quadratic-variable regression
> >
> > > y =3D X*W + e
> >
> > > size(X) =3D [ m n]
> >
> > > n =3D n0*(n0+3)/2
> >
> > > and X contains the columns of X0,
> > > their squares, and their cross products.
> >
> > > Currently, X is created by the following double loop
> >
> > > X =3D X0;
> > > for i =3D 1:n0-1
> > > =A0 =A0 for j=3D i+1:n0
> > > =A0 =A0 =A0 =A0 X =3D [ X X0(:,i).*X0(:,j)];
> > > =A0 =A0 end
> > > end
> > > X =3D [X X0.^2];
> >
> > > Example:
> >
> > > X0 =3D [ =A01 2 3; 4 5 6]
> >
> > > X =3D
> > > =A0 =A0 =A01 =A0 =A0 2 =A0 =A0 3 =A0 =A0 2 =A0 =A0 3
=A0 =A0 6 =A0 =A0 1=
> =A0 =A0 4 =A0 =A0 9
> > > =A0 =A0 =A04 =A0 =A0 5 =A0 =A0 6 =A0 =A020 =A0 =A024
=A0 =A030 =A0 =A016=
> =A0 =A025 =A0 =A036
> >
> > > Can the double loop be vectorized?
> >
> > > TIA,
> >
> > > Greg
> >
> > Hi,
> >
> > It seems that you have messed up the formulation of the
> > quadratic least square problem. Isnt it(?):
> >
> > yi=3Dc0+c1*xi+c2*xi.^2
>
> No.
>
> There are multiple variables ( e.g., n0 =3D 3 ) in X0,
yielding
> n0 original variable columns, n0 variable squared columns
> and n0*(n0-1)/2 crossproduct columns in X.
>
> When n0 =3D 3, n =3D 9 and, for size(y) =3D [m 1]
> size(W) =3D [n+1 1] =3D [10 1]. For measurement i
>
> yhati =3D W1*1 + W2*x1i + W3*x2i + W4*x3i +
> W5*x1i*x2i + W6 *x1i*x3i + W7*x2i*x3i +
> W8*x1i^2 + W9*x2i^2 + W10*x3i^2
>
> where
>
> W =3D [ones(m,1) X]\y;
>
> Hope this helps.
>
> Greg
>
Ah, ok, but this is not much more complicated than I
mentioned before. So your real problem is for example to fit
a rotated and displaced 3D-surface of an ellipsoid to data
points (xi,yi,zi), right? Lets assume X=[x(:) y(:) z(:)];

%--- generate some surface to fit and add noice
clear;clf;
[x,y,z]=ellipsoid(-1/8,-1/32,-1/50,0.5,0.25,0.2,20);
noice=0.01;
X=[x(:) y(:) z(:)]+noice*randn(length(x(:)),3);
figure(1), clf;surf(x,y,z);
hold on;plot3(X(:,1),X(:,2),X(:,3),'.');axis equal;

f=ones(size(x(:))); % specific for this case,
% see equ in: help ellipsoid, i.e., ones are not included.
A=[X(:,1) X(:,2) X(:,3) X(:,1).*X(:,2) X(:,1).*X(:,3) ...
   X(:,2).*X(:,3) X(:,1).^2 X(:,2).^2 X(:,3).^2];
W=A\f

if you set noice=0.0 you get all cross terms equal to zero.

Greg Heath <heath@alumni.brown.edu> wrote in message
<cddc4bd3-1f78-4fad-90d2-a7cde7ec5189@e39g2000hsf.googlegroups.com>...

>
> Can the double loop be vectorized?
>

Yes, for example like this:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X0 = [ 1 2 3; 4 5 6]
[m n0]=size(X0);

[I J]=meshgrid(1:n0,1:n0);

X0_IJ = reshape(X0(:,I).*X0(:,J), m, n0*n0);

itri=[find(J>I) sub2ind(size(I),1:n0,1:n0)'];
I=I(itri); J=J(itri);
k=sub2ind([n0 n0],I,J);

% Bruno
X = [X0 X0_IJ(:,k)]

Slighly more compact

X0 = [ 1 2 3; 4 5 6]
[m n0]=size(X0)

[I J]=meshgrid(1:n0,1:n0);

X0_IJ = reshape(X0(:,I).*X0(:,J), m, n0*n0);

% There might be a better way than using 'find'
itri=[find(J>I); sub2ind(size(I),1:n0,1:n0)'];
I=I(itri); J=J(itri);

X = [X0 X0_IJ(:,sub2ind([n0 n0],I,J))]

% Bruno

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<g4npnp$rsb$1@fred.mathworks.com>...
> Slighly more compact
>
> X0 = [ 1 2 3; 4 5 6]
> [m n0]=size(X0)
>
> [I J]=meshgrid(1:n0,1:n0);
>
> X0_IJ = reshape(X0(:,I).*X0(:,J), m, n0*n0);
>
> % There might be a better way than using 'find'
> itri=[find(J>I); sub2ind(size(I),1:n0,1:n0)'];
> I=I(itri); J=J(itri);
>
> X = [X0 X0_IJ(:,sub2ind([n0 n0],I,J))]
>
> % Bruno

Ok, maby this is a little bit shorter
%--- Prepare some data ---
n=4; % dimensions
Nd=7; % # of data
X=[ones(Nd,1) rand(Nd,n)];
f=rand(Nd,1);

%solve multi-dim second order least square fit
[I,J]=meshgrid(1:n+1);
X=X(:,I(find(I>=J))).*X(:,J(find(I>=J)));
W=X\f;

% present solution
% term: x_i*x_j,Wij (x_0=1).
[int2str(I(find(I>=J))-1) ...
int2str(J(find(I>=J))-1) num2str(W,4)]

On Jul 5, 6:05 am, "Per Sundqvist" <sun...@fy.chalmers.se> wrote:
> Greg Heath <he...@alumni.brown.edu> wrote in message
>
> <6df4a42e-3d83-4dd7-9441-c371361da...@c65g2000hsa.googlegroups.com>...> On Jul 2, 2:16=A0pm, "Per Sundqvist"
> <sun...@fy.chalmers.se> wrote:
> > > Greg Heath <he...@alumni.brown.edu> wrote in message
>
> <cddc4bd3-1f78-4fad-90d2-a7cde7ec5...@e39g2000hsf.googlegroups.com>...
>
> > > > I have an original predictor matrix X0 with
>
> > > > size(X0) =3D [ m n0]
>
> > > > I'd like to perform the
>
> > > > linear-coefficient/quadratic-variable regression
>
> > > > y =3D X*W + e
>
> > > > size(X) =3D [ m n]
>
> > > > n =3D n0*(n0+3)/2
>
> > > > and X contains the columns of X0,
> > > > their squares, and their cross products.
>
> > > > Currently, X is created by the following double loop
>
> > > > X =3D X0;
> > > > for i =3D 1:n0-1
> > > > =A0 =A0 for j=3D i+1:n0
> > > > =A0 =A0 =A0 =A0 X =3D [ X X0(:,i).*X0(:,j)];
> > > > =A0 =A0 end
> > > > end
> > > > X =3D [X X0.^2];
>
> > > > Example:
>
> > > > X0 =3D [ =A01 2 3; 4 5 6]
>
> > > > X =3D
> > > > =A0 =A0 =A01 =A0 =A0 2 =A0 =A0 3 =A0 =A0 2 =A0 =A0 3
>
> =A0 =A0 6 =A0 =A0 1=> =A0 =A0 4 =A0 =A0 9
> > > > =A0 =A0 =A04 =A0 =A0 5 =A0 =A0 6 =A0 =A020 =A0 =A024
>
> =A0 =A030 =A0 =A016=> =A0 =A025 =A0 =A036
>
> > > > Can the double loop be vectorized?
>
> > > > TIA,
>
> > > > Greg
>
> > > Hi,
>
> > > It seems that you have messed up the formulation of the
> > > quadratic least square problem. Isnt it(?):
>
> > > yi=3Dc0+c1*xi+c2*xi.^2
>
> > No.
>
> > There are multiple variables ( e.g., n0 =3D 3 ) in X0,
> yielding
> > n0 original variable columns, n0 variable squared columns
> > and n0*(n0-1)/2 crossproduct columns in X.
>
> > When n0 =3D 3, n =3D 9 and, for size(y) =3D [m 1]
> > size(W) =3D [n+1 1] =3D [10 1]. For measurement i
>
> > yhati =3D W1*1 + W2*x1i + W3*x2i + W4*x3i +
> > W5*x1i*x2i + W6 *x1i*x3i + W7*x2i*x3i +
> > W8*x1i^2 + W9*x2i^2 + W10*x3i^2
>
> > where
>
> > W =3D [ones(m,1) X]\y;
>
> > Hope this helps.
>
> > Greg
>
> Ah, ok, but this is not much more complicated than I
> mentioned before. So your real problem is for example to fit
> a rotated and displaced 3D-surface of an ellipsoid to data
> points (xi,yi,zi), right?

My real problem is multiple regression with m = 1030
observations of n0 = 8 predictors:

size(X0) = [1030 8]
size(X) = [1030 44].

The simple n0 = 3 example was created to clarify my
question.


> Lets assume X=[x(:) y(:) z(:)];
>
> %--- generate some surface to fit and add noice
> clear;clf;
> [x,y,z]=ellipsoid(-1/8,-1/32,-1/50,0.5,0.25,0.2,20);
> noice=0.01;
> X=[x(:) y(:) z(:)]+noice*randn(length(x(:)),3);
> figure(1), clf;surf(x,y,z);
> hold on;plot3(X(:,1),X(:,2),X(:,3),'.');axis equal;
>
> f=ones(size(x(:))); % specific for this case,
> % see equ in: help ellipsoid, i.e., ones are not included.
> A=[X(:,1) X(:,2) X(:,3) X(:,1).*X(:,2) X(:,1).*X(:,3) ...
> X(:,2).*X(:,3) X(:,1).^2 X(:,2).^2 X(:,3).^2];
> W=A\f

That's no help.

My problem is to vectorize the construction of A for
arbitrary n0.

Thanks anyway,

Greg

On Jul 1, 3:05=A0pm, Greg Heath <he...@alumni.brown.edu> wrote:
> I have an original predictor matrix X0 with
>
> size(X0) =3D [ m n0]
>
> I'd like to perform the
>
> linear-coefficient/quadratic-variable regression
>
> y =3D X*W + e
>
> size(X) =3D [ m n]
>
> n =3D n0*(n0+3)/2
>
> and X contains the columns of X0,
> their squares, and their cross products.
>
> Currently, X is created by the following double loop
>
> X =3D X0;
> for i =3D 1:n0-1
> =A0 =A0 for j=3D i+1:n0
> =A0 =A0 =A0 =A0 X =3D [ X X0(:,i).*X0(:,j)];
> =A0 =A0 end
> end
> X =3D [X X0.^2];
>
> Example:
>
> X0 =3D [ =A01 2 3; 4 5 6]
>
> X =3D
> =A0 =A0 =A01 =A0 =A0 2 =A0 =A0 3 =A0 =A0 2 =A0 =A0 3 =A0 =A0 6 =A0 =A0 1 =
=A0 =A0 4 =A0 =A0 9
> =A0 =A0 =A04 =A0 =A0 5 =A0 =A0 6 =A0 =A020 =A0 =A024 =A0 =A030 =A0 =A016 =
=A0 =A025 =A0 =A036
>
> Can the double loop be vectorized?

i was able to reduce it to one loop. However, both
tries took longer than the double loop.

X0 =3D [1 2 3; 4 5 6]
[m n0] =3D size(X0) % [ 1030 8]
tic
for k =3D 1:1e5
    X =3D X0;
    for i =3D 1:n0-1
        for j=3D i+1:n0
            X =3D [ X X0(:,i).*X0(:,j)];
        end
    end
    X =3D [ X X0.^2 ];
end
toc % 2.0630 seconds
tic
for k =3D 1:1e5
    X =3D X0;
    for i =3D 1:n0-1
        X =3D [ X repmat(X0(:,i),1,n0-i).*X0(:,i+1:n0)];
    end
    X =3D [ X X0.^2 ];
end
toc % 14.2810 (factor 6.9224)
tic
for k =3D 1:1e5
    X =3D X0;
    for i =3D 1:n0
        X =3D [ X repmat(X0(:,i),1,n0+1-i).*X0(:,i:n0)];
    end
end
toc % 19.6720 (factor 9.5356)

Greg

Greg Heath <heath@alumni.brown.edu> wrote in message
<6c5247e1-1062-4c25-ab34-12e9ca100527@e39g2000hsf.googlegroups.com>...
> On Jul 1, 3:05=A0pm, Greg Heath <he...@alumni.brown.edu>
wrote:
> > I have an original predictor matrix X0 with
> >
> > size(X0) =3D [ m n0]
> >
> > I'd like to perform the
> >
> > linear-coefficient/quadratic-variable regression
> >
> > y =3D X*W + e
> >
> > size(X) =3D [ m n]
> >
> > n =3D n0*(n0+3)/2
> >
> > and X contains the columns of X0,
> > their squares, and their cross products.
> >
> > Currently, X is created by the following double loop
> >
> > X =3D X0;
> > for i =3D 1:n0-1
> > =A0 =A0 for j=3D i+1:n0
> > =A0 =A0 =A0 =A0 X =3D [ X X0(:,i).*X0(:,j)];
> > =A0 =A0 end
> > end
> > X =3D [X X0.^2];
> >
> > Example:
> >
> > X0 =3D [ =A01 2 3; 4 5 6]
> >
> > X =3D
> > =A0 =A0 =A01 =A0 =A0 2 =A0 =A0 3 =A0 =A0 2 =A0 =A0 3 =A0
=A0 6 =A0 =A0 1 =
> =A0 =A0 4 =A0 =A0 9
> > =A0 =A0 =A04 =A0 =A0 5 =A0 =A0 6 =A0 =A020 =A0 =A024 =A0
=A030 =A0 =A016 =
> =A0 =A025 =A0 =A036
> >
> > Can the double loop be vectorized?
>
> i was able to reduce it to one loop. However, both
> tries took longer than the double loop.
>
> X0 =3D [1 2 3; 4 5 6]
> [m n0] =3D size(X0) % [ 1030 8]
> tic
> for k =3D 1:1e5
> X =3D X0;
> for i =3D 1:n0-1
> for j=3D i+1:n0
> X =3D [ X X0(:,i).*X0(:,j)];
> end
> end
> X =3D [ X X0.^2 ];
> end
> toc % 2.0630 seconds
> tic
> for k =3D 1:1e5
> X =3D X0;
> for i =3D 1:n0-1
> X =3D [ X repmat(X0(:,i),1,n0-i).*X0(:,i+1:n0)];
> end
> X =3D [ X X0.^2 ];
> end
> toc % 14.2810 (factor 6.9224)
> tic
> for k =3D 1:1e5
> X =3D X0;
> for i =3D 1:n0
> X =3D [ X repmat(X0(:,i),1,n0+1-i).*X0(:,i:n0)];
> end
> end
> toc % 19.6720 (factor 9.5356)
>
> Greg

hi, look above message, this X is "A", all is there:
X=X(:,I(find(I>=J))).*X(:,J(find(I>=J)));

Greg Heath <heath@alumni.brown.edu> wrote in message
<6c5247e1-1062-4c25-ab34-12e9ca100527@e39g2000hsf.googlegroups.com>...
>
>
> i was able to reduce it to one loop. However, both
> tries took longer than the double loop.
>


Not sure why you did not include the vectorized code I
proposed in chrono.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X0 =[1 2 3; 4 5 6]
[m n0] =size(X0) % [ 1030 8]
tic
for k =1:1e5
    X =X0;
    for i = 1:n0-1
        for j= i+1:n0
            X = [ X X0(:,i).*X0(:,j)];
        end
    end
    X =[ X X0.^2 ];
end
toc % 4.402068 seconds, Bruno's PC

% I assume these preparation steps can be done before hand
[I J]=meshgrid(1:n0,1:n0);
%
% Or simply itri=find(J>=I)
% if column permutation can be tolerated
%
itri=[find(J>I); sub2ind(size(I),1:n0,1:n0)'];
ij=sub2ind([n0 n0],I(itri),J(itri));

tic
for k =1:1e5
    X0_IJ = reshape(X0(:,I).*X0(:,J), m, n0*n0);
    X = [X0 X0_IJ(:,ij)];
end
toc % 2.567997 seconds, , Bruno's PC

% Bruno

On Jul 5, 8:34 am, "Bruno Luong" <b.lu...@fogale.fr> wrote:
> Slighly more compact
>
> X0 = [ 1 2 3; 4 5 6]
> [m n0]=size(X0)
>
> [I J]=meshgrid(1:n0,1:n0);
>
> X0_IJ = reshape(X0(:,I).*X0(:,J), m, n0*n0);
>
> % There might be a better way than using 'find'
> itri=[find(J>I); sub2ind(size(I),1:n0,1:n0)'];
> I=I(itri); J=J(itri);
>
> X = [X0 X0_IJ(:,sub2ind([n0 n0],I,J))]

This is 18 times slower than using the double loop.

Greg

On Jul 5, 3:34=A0pm, "Per Sundqvist" <sun...@fy.chalmers.se> wrote:
--------SNIP

> %solve multi-dim second order least square fit
> [I,J]=3Dmeshgrid(1:n+1);
> X=3DX(:,I(find(I>=3DJ))).*X(:,J(find(I>=3DJ)));

This is 2.5 times slower than using the double loop.

Greg

On Jul 5, 6:27 pm, "Bruno Luong" <b.lu...@fogale.fr> wrote:
> Greg Heath<he...@alumni.brown.edu> wrote in message
>
> <6c5247e1-1062-4c25-ab34-12e9ca100...@e39g2000hsf.googlegroups.com>...
>
> > i was able to reduce it to one loop. However, both
> > tries took longer than the double loop.
>
> Not sure why you did not include the vectorized code I
> proposed in chrono.

Wha?

> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> X0 =[1 2 3; 4 5 6]
> [m n0] =size(X0) % [ 1030 8]
> tic
> for k =1:1e5
> X =X0;
> for i = 1:n0-1
> for j= i+1:n0
> X = [ X X0(:,i).*X0(:,j)];
> end
> end
> X =[ X X0.^2 ];
> end
> toc % 4.402068 seconds, Bruno's PC

I'm getting from 2 to 4 sec.

> % I assume these preparation steps can be done before hand
> [I J]=meshgrid(1:n0,1:n0);
> %
> % Or simply itri=find(J>=I)
> % if column permutation can be tolerated
> %
> itri=[find(J>I); sub2ind(size(I),1:n0,1:n0)'];
> ij=sub2ind([n0 n0],I(itri),J(itri));
>
> tic
> for k =1:1e5
> X0_IJ = reshape(X0(:,I).*X0(:,J), m, n0*n0);
> X = [X0 X0_IJ(:,ij)];
> end
> toc % 2.567997 seconds, , Bruno's PC

I only have one X0 so I need to include the preparation

tic
for k =1:1e5
     [I J] = meshgrid(1:n0,1:n0);
     itri = [find(J>I); sub2ind(size(I),1:n0,1:n0)'];
     ij = sub2ind([n0 n0],I(itri),J(itri));
     X0_IJ = reshape(X0(:,I).*X0(:,J), m, n0*n0);
     X = [X0 X0_IJ(:,ij)];
 end
 toc % ~ 36 seconds, Greg's PC

and

tic
[I J] = meshgrid(1:n0,1:n0);
itri = [find(J>I); sub2ind(size(I),1:n0,1:n0)'];
ij = sub2ind([n0 n0],I(itri),J(itri));

for k =1:1e5
     X0_IJ = reshape(X0(:,I).*X0(:,J), m, n0*n0);
     X = [X0 X0_IJ(:,ij)];
 end
 toc % ~ 16 seconds, Greg's PC

Not sure why our results differ.

Greg

Greg Heath <heath@alumni.brown.edu> wrote in message
<720bde1b-5f7b-409a-

> I'm getting from 2 to 4 sec.
>

It seems like a main part time is not consumed for generated
triangular indices i and j. The time to compute the vector X
is minimal.

X0 =[1 2 3; 4 5 6]
[m n0] =size(X0) % [ 1030 8]

% Straight forward double-loop
tic
for k =1:1e5
    X =X0;
    for i = 1:n0-1
        for j= i+1:n0
            X = [ X X0(:,i).*X0(:,j)];
        end
    end
    X =[ X X0.^2 ];
end
toc % 4.241084 seconds, Bruno's PC

% Compute the indices i, j by loop
tic
for k =1:1e5
    ii = [];
    jj = [];
    for i = 1:n0-1
        ii = [ii i+zeros(1,n0-i)];
        jj = [jj i+1:n0];
    end
end
toc % 4.248801 seconds, Bruno's PC, almost like above loop

% Compute the indices i, j arrayfun
ifun=@(i) i+zeros(1,n0-i);
jfun=@(i) (i+1:n0);
tic
for k =1:1e5
    ci = arrayfun(ifun,(1:n0),'UniformOutput',false);
    cj = arrayfun(jfun,(1:n0),'UniformOutput',false);
    ii = [ci{:} 1:n0];
    jj = [cj{:} 1:n0];
end
toc % 20.084803 seconds, Bruno's PC

% Compute X array from indices
tic
for k =1:1e5
    X = [X0 X0(:,ii).*X0(:,jj)];
end
toc % 1.850834 seconds, Bruno's PC

% Bruno

Greg Heath <heath@alumni.brown.edu> wrote in message
<cddc4bd3-1f78-4fad-90d2-
a7cde7ec5189@e39g2000hsf.googlegroups.com>...
> I have an original predictor matrix X0 with
>
> size(X0) = [ m n0]
>
> I'd like to perform the
>
> linear-coefficient/quadratic-variable regression
>
> y = X*W + e
>
> size(X) = [ m n]
>
> n = n0*(n0+3)/2
>
> and X contains the columns of X0,
> their squares, and their cross products.
>
> Currently, X is created by the following double loop
>
> X = X0;
> for i = 1:n0-1
> for j= i+1:n0
> X = [ X X0(:,i).*X0(:,j)];
> end
> end
> X = [X X0.^2];
>
> Example:
>
> X0 = [ 1 2 3; 4 5 6]
>
> X =
> 1 2 3 2 3 6 1 4 9
> 4 5 6 20 24 30 16 25 36
>
> Can the double loop be vectorized?
>
> TIA,
>
> Greg

Greg,

Please let me know if this is what you need (the data
matrix here is for easy check of the computation)

% define data matrix X0
X0 = repmat(1:8, [1030, 1]);

% get data size
[m, n0] = size(X0);

% define column index for computation of cross products
% X0(:, [1, 1, ..., 1, 2, 2, ..., 2, ..., n0-2, n0-2,
n0-1])
% .* X0(:, [2, 3, ..., n0, 3, 4, ..., n0, ..., n0-1,
n0, n0])
%
% idea : first factor column index is
% cumsum([1, 0, ..., 0, 1, 0, ..., 0, ..., 1, 0, 1])
where the 1's appear
% at 1, (n0-1)+1, (n0-1)+(n0-2)+1, (n0-1)+(n0-2)+(n0-3)
+1, ...,
% (n0-1)+...+(n0-(n0-2))+1
%
% second factor column index is
% cumsum([2, 1, ..., 1, 3-n0, 1, ..., 1, 4-n0, 1, ...,
n0-n0]) where the
% changes appear at the same places as the 1's in the
first factor column
% index, and the changes are [2, (3:n0)-n0]
%
% compute index of 1's
changeIdx = [1, cumsum(n0-1:-1:2)+1];
% set 1's for first factor column index
ff1(changeIdx) = 1;
% set changes for second factor column index
sfc = ones(1, changeIdx(end));
sfc(changeIdx) = [2, (3:n0)-n0];

% compute matrix of linear system
X = [ones(m, 1), X0, X0(:, cumsum(ff1)).* X0(:, cumsum
(sfc)), X0.^2];

Regards,
Pierre

On Jul 1, 3:05=A0pm, Greg Heath <he...@alumni.brown.edu> wrote:
> I have an original predictor matrix X0 with
>
> size(X0) =3D [ m n0]
>
> I'd like to perform the
>
> linear-coefficient/quadratic-variable regression
>
> y =3D X*W + e
>
> size(X) =3D [ m n]
>
> n =3D n0*(n0+3)/2
>
> and X contains the columns of X0,
> their squares, and their cross products.
>
> Currently, X is created by the following double loop
>
> X =3D X0;
> for i =3D 1:n0-1
> =A0 =A0 for j=3D i+1:n0
> =A0 =A0 =A0 =A0 X =3D [ X X0(:,i).*X0(:,j)];
> =A0 =A0 end
> end
> X =3D [X X0.^2];
>
> Example:
>
> X0 =3D [ =A01 2 3; 4 5 6]
>
> X =3D
> =A0 =A0 =A01 =A0 =A0 2 =A0 =A0 3 =A0 =A0 2 =A0 =A0 3 =A0 =A0 6 =A0 =A0 1 =
=A0 =A0 4 =A0 =A0 9
> =A0 =A0 =A04 =A0 =A0 5 =A0 =A0 6 =A0 =A020 =A0 =A024 =A0 =A030 =A0 =A016 =
=A0 =A025 =A0 =A036
>
> Can the double loop be vectorized?

Many thanks to Tom, Bruno, Per and Pierre.
After some minor modifications, I ran the 4 algorithms
below for my simple 2 X 3 illustrative example and
Pierre's more realistic 1030 X 8 example (the size of
my current regression problem).

Notice how the rankings of the times change for the
two cases. I have noticed that, for the simple example,
the double loop times ranged between 1.8 and 4.2
sec. I can only attribute the spread to whatever
the computer was doing in the backgound without
my permission.

I wouldn't be surprised if the rankings change
when others run the script below.

The following results are for just two runs

run =3D1: 1e5 repetitions for size(X0) =3D [2 3]

 time(sec) Nerd(s)
 2.2970 Greg
 8.1560 Tom
 5.6100 Bruno/Per
 4.1720 Pierre

2. 1e3 repetitions for size(X0) =3D [1030 1]

 time(sec) Nerd(s)
 6.2040 Greg
1.2960 Tom
1.3750 Bruno/Per
1.2040 Pierre % ....The winner by a nose!


%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all, clc

run =3D 2
run =3D 1
if run =3D=3D 1
    X0 =3D [1 2 3; 4 5 6] % [2 3]
    reps =3D 1e5
else
    X0 =3D repmat(1:8,[1030 1]); % [ 1030 8]
    reps =3D 1e3
end
[m n0] =3D size(X0)

tic % Greg
for k =3D 1:reps
    X =3D X0;
    for i =3D 1:n0-1
        for j=3D i+1:n0
            X =3D [ X X0(:,i).*X0(:,j)];
        end
    end
    X =3D [ ones(m,1) X X0.^2 ];
end
t1 =3D toc % 2.3
if run =3D=3D 1, X =3D X, end

% X =3D
% 1 1 2 3 2 3 6 1 4 9
% 1 4 5 6 20 24 30 16 25 36


tic % Tom
for k =3D 1:reps
    col =3D repmat(1:n0,n0,1);
    row =3D col';
    upper =3D ~tril(ones(n0)); % get only (i,j) with j > i
    I =3D row(upper);
    J =3D col(upper);
    X =3D [ones(m,1), X0, X0(:,I).*X0(:,J) X0.^2];
end
t2 =3D toc % ~ 8.2 sec
if run =3D=3D 1, X =3D X, end

tic % Brun0 & Per
for k =3D 1:reps
    [J,I]=3Dmeshgrid(1:n0); % Bruno
    K =3D J>=3DI; % Per
    X =3D [ones(m,1) X0 X0(:,I(K)).*X0(:,J(K))];
end
t3 =3D toc % ~ 5.6 sec
if run =3D=3D 1, X =3D X, end

tic % Pierre
for k =3D 1:reps
  changeIdx =3D [1, cumsum(n0-1:-1:2)+1];
  ff1(changeIdx) =3D 1;
  sfc =3D ones(1, changeIdx(end));
  sfc(changeIdx) =3D [2, (3:n0)-n0];
  X =3D [ones(m,1),X0,X0(:,cumsum(ff1)).*X0(:,cumsum(sfc)),X0.^2];
end
t4=3Dtoc % ~ 4.2 sec
if run =3D=3D 1, X =3D X, end

times =3D [t1 t2 t3 t4]'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Hope this helps.

Greg