mingyan08 / pd3o Goto Github PK

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   min    1/2 x' * A * x + delta(A * x)   (1)
% Where delta(A * x) is Indicator function for \in [a, b]
%
% f(x) = 0.5 * x' * A * x
% g(x) = 0 -> Prox_G(y) = y
% h(A * x) = delta(A x) -> Prox_H(y) = max(min(y, b), a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear
clc;
addpath('fcns','data','output') 
numRows = 5;
numCols = numRows; %<! PSD Matrix
% ---------------------- Generate random data ----------------------
randNum     = 1;
%randn('state', (randNum - 1) * 3 + 1);  % old version
rng((randNum - 1) * 3 + 1, 'v5normal');
A           = rand(numRows, numCols);
A           = A' * A;

% Indicator boundaries: A * x \in [valA, valB]
valA = 0.25;
valB = 1.25;

%----------------------- Set parameters for the problem ------------------------
beta        = 1/normest(A*A');
obj         = PrimalDual;  % using the clas PrimalDual

%% Define all its functions:
% Create handles to functions F, G, H, A, and H*:
obj.myF = @(x) 0.5 * x' * A * x;
obj.myG = @(x) 0;
obj.myH = @(y) 1e6 * any((y > valB) | (y < valA)); %<! Indicator function
obj.myHA = @(y) 0; %<! No use
obj.myA = @(x) A * x;
% Create handles to adjoint/gradient/prox:
obj.myGradF = @(x) A * x;
obj.myProxG = @(x, t) x;
obj.myProxH = @(y, t) min(max(y, valA), valB);
obj.myAdjA  = @(y) A * y; %<! Symmetric

%% Define the parameters for the algorithm
obj.gamma   = 0.5 * beta; % two parameters for the primal-dual algorithms
obj.lambda  = 3/16;       % we will choose different gammas in this example

iter    = 4000;           % the number of iterations

obj.input.x     = A \ (((valA + valB) / 2) * ones(numCols, 1)); %<! Valid initialization
obj.input.s     = zeros(numCols, 1);
obj.input.iter  = iter;


%% Run the primal-dual codes
% PD3O %%%
j   = 1;
tic
[x_PD3O, s_PD3O, E_PD3O, out_PD3O]  = obj.minimize('PD3O', []);
time(j) = toc;

x_PD3O