Revision b5ae1f9b050367711308d270ad23e9d645bee29d authored by Lorenzo Stella on 25 January 2017, 18:59:34 UTC, committed by Lorenzo Stella on 25 January 2017, 18:59:34 UTC
1 parent a2ffa26
forbes_hankel.m
% FORBES_HANKEL
%
% FORBES_HANKEL(U, bx, r, mu, opt) solves
%
% minimize (1/2)||x-b||^2 + mu*sum(svd(H(x)U)
%
% or equivalently
%
% minimize (1/2)||x-b||^2 + mu*sum(svd(z))
% subject to H(x)U - z = 0
%
% Here, H(x) is a block Hankel matrix formed
% from x, and b and U come from measurements:
%
% U is a matrix such that U'U=I
% x is p by N+1 matrix
% U is N-r+1 by n matrix
% H(x) is p(r+1) by N-r+1 Hankel matrix
%
function out = forbes_hankel(U, b, r, mu, tol, freq, opt, out1)
p = size(b,1);
N = size(b,2)-1;
m = (r+1)*p;
n = size(U,2);
f = quadLoss(1, b);
g = nuclearNorm(m, n, mu);
A = opHankel(N, r, p, U);
if nargin < 8 || isempty(out1)
y0 = zeros(m,n);
else
y0 = out1.y;
end
opt.Lf = min(r+1,N-r+1); % squared norm of A
opt.term = @(prob, it, gam, cache_0, cache_x, ops) terminate_hankel(it, freq, cache_x, b, mu, tol);
out = forbes(f, g, y0, [], {A, -1, 0}, opt);
G = -A*out.x1;
nmEY = norm(out.x1-b,'fro')^2/2;
out.pobj = nmEY + mu*sum(svd(G));
end
function flag = terminate_hankel(it, freq, cache, b, mu, tol)
persistent pobj_best dobj_best;
if it == 1
pobj_best = +inf;
dobj_best = -inf;
end
flag = false;
if it <= 1 || mod(it, freq) ~= 0
return;
end
x = cache.gradf1res1x;
G = cache.gradf1x;
nmEY = norm(x-b,'fro')^2/2;
pobj = nmEY + mu*sum(svd(G));
if pobj < pobj_best
pobj_best = pobj;
end
dobj = -cache.FBE;
if dobj > dobj_best
dobj_best = dobj;
end
gap_rel = abs(pobj_best-dobj_best)/max(abs(dobj_best),1);
if gap_rel < tol
flag = true;
end
end
Computing file changes ...