Nine algorithms - sondette.com

I got a request to send the code of the nine algorithms that were compared in “A comparison of nine PLS1 algorithms” and published in the Journal of Chemometrics 2009; 23: 518–529 and led to the Kowalski price 2010 for best theoretical paper. Then I thought it would be even better if they would be easily available for everyone, and that it would be best to publish them here on this blog.

The conclusion of the paper still holds but it could be added that not only deflation will give the best numerical precision. Using modified Gram-Schmidt orthogonalizations on the PLS loading and score vectors will also lead to numerically stable results. On the other hand, when thinking about it, the modified Gram-Schmidt orthogonalization is a kind of deflation too.

One more thing that I have noticed after publishing that could be worth a comment is that the implementations of Krylov PLS and PLSPLS could be as fast as the fastest algorithms if a the mutliplication of X would be stored in a temporary variable, but anyway they will never belong to the top class of numerical precision, so perhaps it does not matter much.

Anyway, here we go with the original code that was published:

NIPALS by Wold

function [b,p,q,w,t] = nipals_pls1(x,y,A);
for a = 1:A,
v = x’*y;
w(:,a) = v/sqrt(v’*v);
t(:,a) = x*w(:,a);
tt = t(:,a)’*t(:,a);
p(:,a) = x’*t(:,a)/tt;
x = x‑t(:,a)*p(:,a)’;
q(a,1) = t(:,a)’*y/tt;
end
b = cumsum(w*(inv(p’*w)*diag(q’)),2);

Non-orthogonalized scores algorithm by Martens

function [b,w,t] = nonorth_scorespls1(x,y,A);
y0 = y;
for a = 1:A
v = x’*y;
w(:,a) = v/sqrt(v’*v);
t(:,a) = x*w(:,a);
qA = inv(t’*t)*t’*y0;
% uncomment for best stability
% qA = t\y0;
b(:,a) = w*qA;
x = x‑t(:,a)*w(:,a)’;
y = y0‑t*qA;
end

Bidiag2 by Golub and Kahan

function [b,w,t,B] = bidiag2(x,y,A);
w(:,1) = x’*y;
w(:,1) = w(:,1)/norm(w(:,1));
t(:,1) = x*w(:,1);
B(1,1) = norm(t(:,1));
t(:,1) = t(:,1)/B(1,1);
for a = 2:A,
w(:,a) = x’*t(:,a-1)-B(a-1,a-1)*w(:,a-1);
B(a-1,a) = norm(w(:,a));
w(:,a) = w(:,a)/B(a-1,a);
t(:,a) = x*w(:,a)-B(a-1,a)*t(:,a-1);
B(a,a) = norm(t(:,a));
t(:,a) = t(:,a)/B(a,a);
end
invB = inv(B);
q = y’*t;
for a = 1:A,
b(:,a) = w(:,1:a)*(invB(1:a,1:a)*q(1:a)’);
end

SIMPLS by de Jong

function [b,R,T,P,Q] = simpls1(x,y,A)
[tmp,K] = size(x);
V = zeros(K,A);
s = x’*y;
for a = 1:A
r = s;
t = x*r;
tt = norm(t);
t = t/tt;
r = r/tt;
p = x’*t;
q = y’*t;
u = y*q;
v = p;
if a > 1
v = v — V*(V’*p);
u = u — T*(T’*u);
end
v = v/norm(v);
s = s — v*(v’*s);
R(:,a) = r;
T(:,a) = t;
P(:,a) = p;
Q(:,a) = q;
U(:,a) = u;
V(:,a) = v;
end
b = cumsum(R*diag(Q),2);

Improved kernel PLS by Dayal

function [b,W,P,T,Q,R] = impr_kernel1(X,Y,A);
W=[];R=[];P=[];Q=[];T=[];U=[];
XY=X’*Y;
for a=1:A
w=XY;
w=w/sqrt(w’*w);
r=w;
for k=1:a-1,
r=r-(P(:,k)’*w)*R(:,k);
end
t = X*r;
tt = t’*t;
p = X’*t/tt;
q = (r’*XY)’/tt;
XY = XY-(p*q’)*tt;
W = [W w];
P = [P p];
T = [T t];
Q = [Q q];
R = [R r];
end
b = cumsum(R*diag(Q’),2);

PLSF by Manne

function b = plsf(x,y,A);
W(:,1) = x’*y;
nw1 = norm(W(:,1));
W(:,1) = W(:,1)/nw1;
U = x*W(:,1);
Rd(1) = norm(U);
U = U/Rd(1);
Uy(1) = U’*y;
k = 1;
while k > 0,
k = k+1;
W(:,k) = x’*U‑W(:,k-1)*Rd(k-1);
Ru(k-1) = norm(W(:,k));
W(:,k) = W(:,k)/Ru(k-1);
if k > A,
break;
else
U = x*W(:,k)-U*Ru(k-1);
Rd(k) = norm(U);
U = U/Rd(k);
Uy(k) = U’*y;
end
end
W(:,k) = [];
Ru(k-1) = [];
n1 = k-1;
R = zeros(n1,n1);
R = R+diag(Rd)+diag(Ru,1);
Uy = Uy’;
invR = inv(R);
for k = 1:A,
b(:,k) = W(:,1:k)*(invR(1:k,1:k)*Uy(1:k));
end

Direct-scores PLS by Andersson

function [b,U,R,P,q] = directscorespls1(x,y,A)
s = x’*y;
utemp = x*s;
nrm = norm(utemp);
R(:,1) = s/nrm;
U(:,1) = utemp/nrm;
P(:,1) = x’*U(:,1);
q(:,1) = y’*U(:,1);
b(:,1) = R*q’;
for a = 2:A,
s = s‑P(:,a-1)*q(:,a-1)’;
rtemp = s — R*(P’*s);
utemp = x*rtemp;
nrm = norm(utemp);
R(:,a) = rtemp/nrm;
U(:,a) = utemp/nrm;
P(:,a) = x’*U(:,a);
q(:,a) = y’*U(:,a);
b(:,a) = R*q’;
end

Krylov PLS1 by Andersson

function b = krylovpls1(x,y,A)
V = x’*y;
b = V*((x*V)\y);
for a = 2:A
V = [V x’*(x*V(:,a-1))];
b(:,a) = V*((x*V)\y);
end

PLSPLS1 by Andersson

function b = plspls1(x,y,A)
V = x’*y;
b = V*((x*V)\y);
for a = 2:A
V = [b x’*(x*b(:,a-1))];
b(:,a) = V*((x*V)\y);
end

Leave a Comment Cancel Reply