Skip to content
Snippets Groups Projects
Commit 1d162ea8 authored by Duong Hung Pham's avatar Duong Hung Pham
Browse files

Upload New File

parent a5b47b04
Branches
No related tags found
No related merge requests found
SRPCA.m 0 → 100644
function [y, x, errxn] = SRPCA(S, H, lambda, mu, max_iter)
% DRPCA.
%
% Input
% - S is a data matrix (of the size n x m) to be decomposed
% S can also contain NaN's for unobserved values
% - H: Point Spread function
% - lambda - sparse regularization parameter, default = 1/sqrt(max(m*n,p))
% - mu - low-rank regularization parameter, default = 10
% - max_iter - maximum number of iterations, default = 20
% - rho - convergence parameter, default = 1e-6
% - eps - accelerate convergence
% - rho_max - maximum rho
% - tol - reconstruction error tolerance, default = 1e-6
%
% Ouput
%- x is high resolution blood
% - T is tissue
[m, n, p] = size(S);
unobserved = isnan(S);
S(unobserved) = 0;
normX = norm(S(:), 'fro');
if nargin > 1 && any(H(:))
[Mh, Nh] = size(H);
center = round([Mh, Nh] / 2);
H = fft2(circshift(padarray(H, [m - Mh, n - Nh], 'post'), 1 - center));
else
H = ones(m,n);
end
if nargin < 3
lambda = 10/sqrt(m*n);
end
if nargin < 4
mu = 10;
end
if nargin < 5
max_iter = 20;
end
rho = 1e-6;
eps = 1.5;
rho_max = 0.0033;
tol = 1e-6;
%% initial solution
T = zeros(m, n, p);
y = zeros(m, n, p);
x = zeros(m, n, p);
Z = zeros(m, n, p);
N = zeros(m, n, p); % gamma1
W = zeros(m, n, p); % gamma2
Hx = zeros(m, n, p);
errxn= zeros(1,max_iter);
Dt = conj(H);
DD = abs(H).^2;
for iter = (1:max_iter)
xold=x;
y = (S-Hx + rho*(T - W))./(1+rho);
x = ifft2(fft2(ifft2(Dt.*fft2(S-y)) + rho*(Z - N))./(DD + rho));
Hx = ifft2(H.*fft2(x));
Z = So(lambda, x + N);
[T, RankT] = Do(mu, reshape(y + W,m*n,p));
T = reshape(T,m,n,p);
Z1 = y - T;
Z1(unobserved) = 0; % skip missing values
W = W + Z1; %gamma2
Z2 = x - Z;
Z2(unobserved) = 0; % skip missing values
N = N + Z2; %gamma1
rho = min(rho*eps, rho_max);
Z3 = S - T - Hx;
err1 = norm(Z1(:), 'fro') / normX;
err2 = norm(Z2(:), 'fro') / normX;
err3 = norm(Z3(:), 'fro') / normX;
errxn(1,iter) = norm(xold(:)-x(:), 2) / norm(xold(:), 2);
if (iter == 1) || (err1 > tol) || (err2 > tol) || (err3 > tol)
fprintf(1, 'iter: %04d\terr1: %f\terr2: %f\terr3: %f\trank(T): %d\tcard(S): %d\n', ...
iter, err1, err2, err3, RankT, nnz(x(~unobserved)));
end
if (err1 < tol) && (err2 < tol)
break;
end
end
end
function r = So(tau, S)
% shrinkage operator
r = sign(S) .* max(abs(S) - tau, 0);
end
function [r, RankT] = Do(tau, S)
% shrinkage operator for singular values
[U, D, V] = svd(S, 'econ');
r = U*So(tau, D)*V';
RankT = rank(r);
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment