function P = PPLS_sift_gen(X, pub, pri, ncomp, lambda)
% Privacy Partial Least Squares (PPLS)
%
% Input:
% X: nxd matrix, n is the number of samples and d is the dimensionality
% pub: NxD1 matrix, indicating the presence/absence of the D1 public attributes (via binary labels)
% pri: NxD1 matrix, indicating the presence/absence of the D1 private attributes (via binary labels)
% ncomp: the maximum sift dimensionality desired (default is 10)
% lamda: tradeoff parameter between privacy/utility used in the objective function (default is 1)
% Output:
% P: sift (transformation) matrix
% Notes:
% This code is code modified from the SIMPLS subfunction in PLSREGRESS
% http://www.mathworks.com/help/stats/plsregress.html
%
% If you use this code please cite the following paper:
% SensorSift: Balancing Sensor Data Privacy and Utility
% in Automated Face Understanding,
% Miro Enev, Jaeyeon Jung, Liefeng Bo, Xiaofeng Ren, Tadayoshi Kohno
% Annual Computer Security Applications Conference (ACSAC), 2012
% Copyright (c) 2013, Miro Enev, Jaeyeon Jung, Liefeng Bo, Xiaofeng Ren, Tadayoshi Kohno
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
% 1. Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
% 3. All advertising materials mentioning features or use of this software
% must display the following acknowledgement:
% This product includes software developed by the Intel Science and
% Technology Center.
% 4. Neither the name of the Intel Science and Technology Center nor the
% names of its contributors may be used to endorse or promote products
% derived from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY US ''AS IS'' AND ANY
% EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
% ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
if nargin < 4, ncomp = size(pub, 2); end
if nargin < 5, lambda = 10; end
meanX = mean(X,1);
X0 = bsxfun(@minus, X, meanX);
P = simpls(X0,pub,pri,ncomp,lambda);
function P = simpls(X0,pub,pri,ncomp,lambda)
[n,dx] = size(X0);
dy = size(pub,2);
% Preallocate outputs
outClass = superiorfloat(X0,pub);
Xloadings = zeros(dx,ncomp,outClass);
% An orthonormal basis for the span of the X loadings, to make the successive
% deflation X0'*pub simple - each new basis vector can be removed from Cov
% separately.
V = zeros(dx,ncomp);
Cov = X0'*pub;
priCovsquare = (X0'*pri)*pri'*X0;
P = [];
for i = 1:ncomp
% Find unit length ti=X0*ri and ui=pub*ci whose covariance, ri'*X0'*pub*ci, is
% jointly maximized, subject to ti'*tj=0 for j=1:(i-1).
[ri,si,ci] = svd(Cov,'econ'); ri = ri(:,1); ci = ci(:,1); si = si(1);
% find transformation matrix
CovSquare = ( Cov*Cov' - lambda*priCovsquare );
[ri, beta] = eigs( CovSquare, 1, 'LR');
% sort eigenvectors (sift components)
ri = real( ri );
beta = real( diag( beta ) );
[~, order] = sort( beta, 'descend' );
% append h to sift matrix
ri = ri( :, order(1) );
P = [P ri];
ti = X0*ri;
normti = norm(ti);
ti = ti ./ normti; % ti'*ti == 1
Xloadings(:,i) = X0'*ti;
% Update the orthonormal basis with modified Gram Schmidt (more stable),
vi = Xloadings(:,i);
for repeat = 1:2
for j = 1:i-1
vj = V(:,j);
vi = vi - (vi'*vj)*vj;
end
end
vi = vi ./ norm(vi);
V(:,i) = vi;
% Deflate Cov, i.e. project onto the ortho-complement of the X loadings.
% First remove projections along the current basis vector, then remove any
% component along previous basis vectors that's crept in as noise from
% previous deflations.
Cov = Cov - vi*(vi'*Cov);
Vi = V(:,1:i);
Cov = Cov - Vi*(Vi'*Cov);
disp([i ncomp]);
end