Contents

Regularized matrix linear regression

clear;
% reset random seed
s = RandStream('mt19937ar','Seed',2);
RandStream.setGlobalStream(s);

True coefficients for regular (non-array) covariates

p0 = 5;
b0 = ones(p0,1);

2D true signal: 64-by-64 cross

shape = imread('cross.gif');
shape = array_resize(shape,[32,32]); % 32-by-32
b = zeros(2*size(shape));
b((size(b,1)/4):(size(b,1)/4)+size(shape,1)-1, ...
    (size(b,2)/4):(size(b,2)/4)+size(shape,2)-1) = shape;
[p1,p2] = size(b);
display(size(b));
ans =

    64    64

Simulate covariates

n = 500;    % sample size
X = randn(n,p0);   % n-by-p0 regular design matrix
M = tensor(randn(p1,p2,n));  % p1-by-p2-by-n matrix variates
display(size(M));
ans =

    64    64   500

Simulate responses

mu = X*b0 + double(ttt(tensor(b), M, 1:2));
sigma = 1;  % noise level
y = mu + sigma*randn(n,1);

Determine max lambda to start

[~,~,stats] = matrix_sparsereg(X,M,y,inf,'normal');
maxlambda = stats.maxlambda*.95;

Fit nuclear norm regularized linear regression at grid points

gridpts = 10;
lambdas = zeros(1,gridpts);
gs = 2/(1+sqrt(5));
B = cell(1,gridpts);
AIC = zeros(1,gridpts);
BIC = zeros(1,gridpts);
tic;
for i=1:gridpts
    if (i==1)
        B0 = [];
    else
        B0 = B{i-1};    % warm start
    end
    lambda = maxlambda*gs^(i-1);
    lambdas(i) = lambda;
    [beta0,B{i},stats] = matrix_sparsereg(X,M,y,lambda,'normal','B0',B0);
    AIC(i) = stats.AIC;
    BIC(i) = stats.BIC;
end
toc;
Elapsed time is 1.597956 seconds.

Display true signal and snapshots along nuclear norm solution path

figure; hold on;
set(gca,'FontSize',20);

subplot(2,2,1);
imagesc(-b);
colormap(gray);
title('True Signal');
axis equal;
axis tight;

ploti = 1;
for i=[gridpts round(gridpts/2) 1]
    ploti = ploti + 1;
    subplot(2,2,ploti);
    imagesc(-double(B{i}));
    colormap(gray);
    title({['nuclear norm,', ' \lambda=', ...
        num2str(lambdas(i))]; ['BIC=', num2str(BIC(i))]});
    axis equal;
    axis tight;
end

Display AIC/BIC trace plot

figure;
set(gca,'FontSize',20);
semilogx(lambdas, AIC, '-+', lambdas, BIC, '-o');
xlabel('\lambda');
ylabel('BIC');
xlim([min(lambdas) max(lambdas)]);
title('Nuclear norm AIC/BIC');
legend('AIC', 'BIC', 'Location', 'northwest');

Compare to lasso penalized linear regression

Transform matrix variates to vector form

TM = tenmat(tensor(M),3,[1 2]);
Xall = [double(TM) X];

Determine max lambda to start

lambdastart = max(lsq_maxlambda(sum(Xall.^2),-y'*Xall,'enet',1));
maxlambda_lasso = lambdastart*.95;

Optimization at grid points

gridpts = 10;
B_lasso = cell(1,gridpts);
BIC_lasso = zeros(1,gridpts);
lambdas_lasso = zeros(1,gridpts);
penidx = true(size(Xall,2),1);
penidx(numel(b)+1:end) = false; % do not penalize regular covariates

tic;
for i=1:gridpts
    if (i==1)
        x0 = zeros(size(Xall,2),1);
    else
        x0 = beta;
    end
    lambda = maxlambda_lasso*gs^(i-1);
    lambdas_lasso(i) = lambda;
    [beta] = lsq_sparsereg(Xall,y,lambda,'x0',x0,'penidx',penidx);
    B_lasso{i} = reshape(beta(1:numel(b)),p1,p2);
    BIC_lasso(i) = .5*norm(y-Xall*beta)^2+log(n)*nnz(abs(beta)>1e-8);
end
toc;
Elapsed time is 0.433877 seconds.

Display true signal and snapshots along lasso solution path

figure; hold on;
set(gca,'FontSize',20);

subplot(2,2,1);
imagesc(-b);
colormap(gray);
title('True Signal');
axis equal;
axis tight;

ploti = 1;
for i=[gridpts round(gridpts/2) 1]
    ploti = ploti + 1;
    subplot(2,2,ploti);
    imagesc(-B_lasso{i});
    colormap(gray);
    title({['lasso' ', \lambda=', ...
        num2str(lambdas(i))]; ['BIC=', num2str(BIC(i))]});
    axis equal;
    axis tight;
end

Lasso BIC path

figure;
set(gca,'FontSize',20);
semilogx(lambdas_lasso,BIC_lasso,'-o');
title('lasso BIC');
xlabel('\lambda');
ylabel('BIC');
xlim([min(lambdas_lasso),max(lambdas_lasso)]);

Regularized matrix Poisson (log-linear) regression

clear;
% reset random seed
s = RandStream('mt19937ar','Seed',2);
RandStream.setGlobalStream(s);

2D true signal: 64-by-64 cross

shape = imread('cross.gif');
shape = array_resize(shape,[32,32]); % 32-by-32
b = zeros(2*size(shape));
b((size(b,1)/4):(size(b,1)/4)+size(shape,1)-1, ...
    (size(b,2)/4):(size(b,2)/4)+size(shape,2)-1) = shape;
[p1,p2] = size(b);
display(size(b));
ans =

    64    64

True coefficients for regular (non-array) covariates

p0 = 5;
b0 = ones(p0,1);

Simulate covariates

n = 750;    % sample size
X = randn(n,p0);   % n-by-p regular design matrix
M = tensor(randn(p1,p2,n));  % n p1-by-p2 matrix variates
display(size(M));

% Simulate Poisson count responses from the systematic components
mu = X*b0 + double(ttt(tensor(b), M, 1:2));
mu = mu/(max(abs(mu)))*5;   % scale to [-5, 5], to avoid overflow
y = poissrnd(exp(mu));
ans =

    64    64   750

Determine max lambda to start

[~,~,stats] = matrix_sparsereg(X,M,y,inf,'poisson');
maxlambda = stats.maxlambda*.95;

Fit nuclear norm regularized Poisson regression at grid points

gridpts = 10;
lambdas = zeros(1,gridpts);
gs = 2/(1+sqrt(5));
B = cell(1,gridpts);
AIC = zeros(1,gridpts);
BIC = zeros(1,gridpts);
tic;
for i=1:gridpts
    if (i==1)
        B0 = [];
    else
        B0 = B{i-1};
    end
    lambda = maxlambda*gs^(i-1);
    lambdas(i) = lambda;
    [beta0,B{i},stats] = matrix_sparsereg(X,M,y,lambda,'poisson','B0',B0);
    AIC(i) = stats.AIC;
    BIC(i) = stats.BIC;
end
toc
Elapsed time is 3.958110 seconds.

Display true signal and snapshots along nuclear norm solution path

figure; hold on;
set(gca,'FontSize',20);

subplot(2,2,1);
imagesc(-b);
colormap(gray);
title('True Signal');
axis equal;
axis tight;

ploti = 1;
for i=[gridpts round(gridpts/2) 1]
    ploti = ploti + 1;
    subplot(2,2,ploti);
    imagesc(-double(B{i}));
    colormap(gray);
    title({['nuclear norm,', ' \lambda=', ...
        num2str(lambdas(i))]; ['BIC=', num2str(BIC(i))]});
    axis equal;
    axis tight;
end

Display AIC/BIC trace plot

figure;
set(gca,'FontSize',20);
semilogx(lambdas, AIC, '-+', lambdas, BIC, '-o');
xlabel('\lambda');
ylabel('BIC');
xlim([min(lambdas) max(lambdas)]);
legend('AIC', 'BIC', 'Location', 'northwest');

Compare to lasso penalized Poisson (log-linear) regression

Transform matrix variates to vector form

TM = tenmat(tensor(M),3,[1 2]);
Xall = [double(TM) X];

Determine max lambda to start

lambdastart = 0;
for j=1:numel(b)
    lambdastart = max(lambdastart,glm_maxlambda(Xall(:,j),y,'loglinear'));
end
maxlambda_lasso = lambdastart*.95;

Optimization at grid points

gridpts = 10;
B_lasso = cell(1,gridpts);
BIC_lasso = zeros(1,gridpts);
lambdas_lasso = zeros(1,gridpts);
penidx = true(size(Xall,2),1);
penidx(numel(b)+1:end) = false; % do not penalize regular covariates

tic;
for i=1:gridpts
    if (i==1)
        x0 = zeros(size(Xall,2),1);
    else
        x0 = beta;
    end
    lambda = maxlambda_lasso*gs^(i-1);
    lambdas_lasso(i) = lambda;
    [beta] = glm_sparsereg(Xall,y,lambda,'loglinear', ...
        'x0',x0,'penidx',penidx);
    B_lasso{i} = beta(1:numel(b));
    eta = Xall*beta;
    BIC_lasso(i) = - sum(y.*log(eta)-gammaln(y+1)-eta) ...
        + log(n)*nnz(abs(beta)>1e-8);
end
toc;
Elapsed time is 66.000598 seconds.

Display true sinal and snapshots along lasso solution path

figure; hold on;
set(gca,'FontSize',20);

subplot(2,2,1);
imagesc(-b);
colormap(gray);
title('True Signal');
axis equal;
axis tight;

ploti = 1;
for i=[gridpts round(gridpts/2) 1]
    ploti = ploti + 1;
    subplot(2,2,ploti);
    imagesc(-reshape(B_lasso{i},p1,p2));
    colormap(gray);
    title({['lasso' ', \lambda=', ...
        num2str(lambdas(i))]; ['BIC=', num2str(BIC(i))]});
    axis equal;
    axis tight;
end

Lasso BIC path

figure;
set(gca,'FontSize',20);
semilogx(lambdas_lasso,BIC_lasso,'-o');
title('lasso BIC');
xlabel('\lambda');
ylabel('BIC');
xlim([min(lambdas_lasso),max(lambdas_lasso)]);
Warning: Imaginary parts of complex X and/or Y arguments ignored 

Regularized matrix logistic regression

clear;
% reset random seed
s = RandStream('mt19937ar','Seed',2);
RandStream.setGlobalStream(s);

2D true signal: 64-by-64 cross

shape = imread('cross.gif');
shape = array_resize(shape,[32,32]); % 32-by-32
b = zeros(2*size(shape));
b((size(b,1)/4):(size(b,1)/4)+size(shape,1)-1, ...
    (size(b,2)/4):(size(b,2)/4)+size(shape,2)-1) = shape;
[p1,p2] = size(b);
display(size(b));
ans =

    64    64

True coefficients for regular (non-array) covariates

p0 = 5;
b0 = ones(p0,1);

Simulate covariates

n = 1000;    % sample size
X = randn(n,p0);   % n-by-p regular design matrix
M = tensor(randn(p1,p2,n));  % n p1-by-p2 matrix variates
display(size(M));
ans =

          64          64        1000

Simulate binary responses from the systematic components

mu = X*b0 + double(ttt(tensor(b), M, 1:2));
y = binornd(1, 1./(1+exp(-mu)));

Determine max lambda to start

[~,~,stats] = matrix_sparsereg(X,M,y,inf,'binomial');
maxlambda = stats.maxlambda*0.95;

Fit nuclear norm regularized logistic regression at grid points

gridpts = 10;
lambdas = zeros(1,gridpts);
gs = 2/(1+sqrt(5));
B = cell(1,gridpts);
AIC = zeros(1,gridpts);
BIC = zeros(1,gridpts);

tic;
for i=1:gridpts
    if (i==1)
        B0 = [];
    else
        B0 = B{i-1};
    end
    lambda = maxlambda*gs^(i-1);
    lambdas(i) = lambda;
    [beta0,B{i},stats] = matrix_sparsereg(X,M,y,lambda,'binomial','B0',B0);
    AIC(i) = stats.AIC;
    BIC(i) = stats.BIC;
end
toc
Elapsed time is 2.765088 seconds.

Display true signal and and snapshots along nuclear norm solution path

figure; hold on;
set(gca,'FontSize',20);

subplot(2,2,1);
imagesc(-b);
colormap(gray);
title('True Signal');
axis equal;
axis tight;

ploti = 1;
for i=[gridpts round(gridpts/2) 1]
    ploti = ploti + 1;
    subplot(2,2,ploti);
    imagesc(-double(B{i}));
    colormap(gray);
    title({['nuclear norm,', ' \lambda=', ...
        num2str(lambdas(i))]; ['BIC=', num2str(BIC(i))]});
    axis equal;
    axis tight;
end

Display AIC/BIC trace plot

figure;
set(gca,'FontSize',20);
semilogx(lambdas, AIC, '-+', lambdas, BIC, '-o');
xlabel('\lambda');
ylabel('BIC');
xlim([min(lambdas) max(lambdas)]);
title('Nuclear norm AIC/BIC');
legend('AIC', 'BIC', 'Location', 'northwest');

Compare to lasso penalized logistic regression

Transform matrix variates to vector form

TM = tenmat(tensor(M),3,[1 2]);
Xall = [double(TM) X];

Determine max lambda to start

lambdastart = 0;            % find the maximum tuning parameter to start
for j=1:numel(b)
    lambdastart = max(lambdastart,glm_maxlambda(Xall(:,j),y,'logistic'));
end
maxlambda_lasso = lambdastart;

Optimization at grid points

gridpts = 10;
B_lasso = cell(1,gridpts);
BIC_lasso = zeros(1,gridpts);
lambdas_lasso = zeros(1,gridpts);
penidx = true(size(Xall,2),1);
penidx(numel(b)+1:end) = false; % do not penalize regular covariates

tic;
for i=1:gridpts
    if (i==1)
        x0 = zeros(size(Xall,2),1);
    else
        x0 = beta;
    end
    lambda = maxlambda_lasso*gs^(i-1);
    lambdas_lasso(i) = lambda;
    [beta] = glm_sparsereg(Xall,y,lambda,'logistic','x0',x0, ...
        'penidx',penidx);
    B_lasso{i} = beta(1:numel(b));
    eta = Xall*beta;
    BIC_lasso(i) = - sum(y.*eta-log(1+exp(eta))) ...
        + log(n)*nnz(abs(beta)>1e-8);
end
toc;
Elapsed time is 111.274021 seconds.

Display true signal and snapshots along lasso solution path

figure; hold on;
set(gca,'FontSize',20);

subplot(2,2,1);
imagesc(-b);
colormap(gray);
title('True Signal');
axis equal;
axis tight;

ploti = 1;
for i=[gridpts round(gridpts/2) 1]
    ploti = ploti + 1;
    subplot(2,2,ploti);
    imagesc(-reshape(B_lasso{i},p1,p2));
    colormap(gray);
    title({['lasso' ', \lambda=', ...
        num2str(lambdas(i))]; ['BIC=', num2str(BIC(i))]});
    axis equal;
    axis tight;
end

Lasso BIC path

figure;
set(gca,'FontSize',20);
semilogx(lambdas_lasso,BIC_lasso,'-o');
title('lasso BIC');
xlabel('\lambda');
ylabel('BIC');
xlim([min(lambdas_lasso),max(lambdas_lasso)]);