Path: news.mathworks.com!not-for-mail
From: "Per Sundqvist" <sunkan@fy.chalmers.se>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Vectorization for Quadratic Polynomial Regression
Date: Sat, 5 Jul 2008 10:05:03 +0000 (UTC)
Organization: Chalmers Tekniska H&#246;gskola
Lines: 108
Message-ID: <g4nh0f$15o$1@fred.mathworks.com>
References: <cddc4bd3-1f78-4fad-90d2-a7cde7ec5189@e39g2000hsf.googlegroups.com>  <6df4a42e-3d83-4dd7-9441-c371361dab06@c65g2000hsa.googlegroups.com>
Reply-To: "Per Sundqvist" <sunkan@fy.chalmers.se>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1215252303 1208 172.30.248.38 (5 Jul 2008 10:05:03 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sat, 5 Jul 2008 10:05:03 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 266682
Xref: news.mathworks.com comp.soft-sys.matlab:477629


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.