# to unpack, sh this message
PATH=/bin:/usr/bin
echo "Anything free comes with no guarantee!"
====================README==========================
This dir contains the matlab codes for the paper "a reduced hessian
method for constrained optimization" by Aiping Liao 
accepted by Computational Optimization and Applications on 4/2/93.
The main program is "secant.m" 
----Aiping Liao
----ACRI, Cornell, U
----Ithaca, NY 14853
----(607)254-8838
----June 28, 1993
===================END README======================
===================secant.m========================
function [xopt,f,k,cv] = secant(x, kmax)
%
% this program is to solve 
% min f(x)
% c_{i}(x) = 0
% c_{j}(x) >= 0.
%
% The user should provide "fgv.m" which evaluates the value and gradient
% of the objective and "constr.m" which provides the value and Jacobian
% of the constraints.
%
% On input:
%       x: the initial point.
%       kmax: the max number of iteration set by user.
% On return:
%       xopt: the solution.
%       f: the objective value at xopt.
%       k: the total number of iterations that took.
%       cv: the constraint-violation.  
%
[n,lo] = size(x);
gg=zeros(n,1);
g=zeros(n,1);
ct = zeros(4,1);
k=0
[f,g] = fgv(x);
ct(1) = ct(1) + 1
ct(2) = ct(2) + 1;
%
%  construct the active space
%
[vc, DC, m1, m] = constr(x);
ct(3) = ct(3) + 1;
ct(4) = ct(4) + 1;
A = DC(:,1:m1);
ca=vc(1:m1);
ai=[1:m1];
if (m1 < m)
 for i = m1+1:m;
   if (abs(vc(i)) < 0.0001)
 A = [A,DC(:,i)];
 ca = [ca;vc(i)];
 ai = [ai,i];
   end
 end
end
[nnt,t] = size(A);
if (t == 0)|(t >= n)
Z = diag(ones(n,1));
else
[Q,R]=qr(A);
Z = Q(:,[t+1:n]);
end
%
B = Z'*diag(ones(n,1))*Z;
H = Z*B*Z';
%
r = 100*ones(m,1);
d=ones(n,1);
[w, dw] = wmerit(x,d,r);
dw = ones(n,1);
while ( (k < kmax) & (norm(dw) > 10^(-9)))
%while (k < kmax)
k = k+ 1
[f,g]=fgv(x);
[vc, DC, m1, m] = constr(x);
if (t < n)
H = Z*B*Z';
%else
%H = diag(ones(n,1));
end
[d, lambda] = qp(H+0*diag(ones(n,1)), g, -DC', vc, [], [], zeros(n,1), m1, [])
dl0 = dlag(x, lambda);
r = abs(lambda) + 10;
%r = 100*ones(m,1);
%
% line search
ls0=1.0;
sigma = 0.0001;
SC = - A*(A'*A)^(-1);
[vc1, DC1, m1, m] = constr(x+d);
d1 = SC*vc1(ai);
for j = 0:20
    ls = ls0*(0.5^j);
    xtem = x + ls*d + ls*ls*d1;
[wtem, dwtem] = wmerit(xtem,d,r);
    if (wtem <= w + sigma*ls*dwtem)
     break;
    end
end
ct(1) = ct(1) + j + 1;
ct(2) = ct(2) + 1;
ct(3) = ct(3) + j + 1;
ct(4) = ct(4) + 1;
ls
s = ls*Z'*d;
dl1 = dlag(x + ls*Z*Z'*d, lambda);
ct(2) = ct(2) + 1;
ct(4) = ct(4) + 1;
y=Z'*(dl1 - dl0);
x=xtem
f=fgv(x)
w=wtem;
% %%%%%%
if (s'*y > 0.2*s'*B*s)
 theta = 1.0;
else
 theta = (0.8*s'*B*s)/(s'*B*s - s'*y);
end
yy = theta*y + (1-theta)*B*s;
B=bfgs(s, yy, B);
%
% identify the active space
%
Zold = Z;
[vc, DC, m1, m] = constr(x);
A = DC(:,1:m1);
ca = vc(1:m1);
ai = [1:m1];
if (m1 < m)
 for i = m1+1:m;
   if (abs(vc(i)) < 0.0001)|(abs(lambda(i)) > 0.0001)
 A = [A,DC(:,i)];
 ca = [ca;vc(i)];
 ai=[ai,i];
   end
 end
end
[nnt,t] = size(A);
% if A if n by n, do not use update: set Z=I
if (t < n)
  if (t > 0)
[Q,R]=qr(A);
Z = Q(:,[t+1:n]);
  else
Z = diag(ones(n,1));
  end
else
Z = diag(ones(n,1));
end
T = Zold'*Z;
[e1,e2]=size(B);
if (t < n)
B = T'*(B - diag(ones(e1,1)))*T + diag(ones(n-t,1));
else
B = T'*(B - diag(ones(e1,1)))*T + diag(ones(n,1));
end
% end
H=B;
norm(B)
dw = dwtem;
cv = cmerit(x)
end
xopt=x;
=================END secant.m==================
=================bfgs.m========================
function B = bfgs(s, y, B)
sbs = s'*B*s;
ys = y'*s;
if (sbs > 10^(-100)) & (abs(ys) > 10^(-100))
B=B-((B*s)*(B*s)')/(s'*B*s) + (y*y')/(y'*s);
end
end
================END bfgs.m====================
================cmerit.m======================
function cv = cmerit(x)
[n,nx] = size(x);
[vc,DC,m1,m] = constr(x);
cv = 0;
if (m1 > 0)
 for i=1:m1;
  cv=cv + (abs(vc(i)));
 end
end
     if (m1 < m)
for i=m1+1:m;
tt=min([0;vc(i)]);
cv = cv - tt;
end
     end
end
===============END cmerit.m==================
===============wmerit.m======================
function [w, dw] = wmerit(x,d,r)
[n,nx] = size(x);
[m,mr] = size(r);
[f,g] =fgv(x);
[vc,DC,m1,m] = constr(x);
w = f;
dw = g'*d;
if (m1 > 0)
 for i=1:m1;
  w=w + r(i)*(abs(vc(i)));
  dw = dw - r(i)*(abs(vc(i)));
 end
end
      if (m1 < m)
for i=m1+1:m;
tt=min([0;vc(i)]);
w = w - r(i)*tt;
end
for i=m1+1:m;
 if (vc(i) < 0)
    fi = d'*DC(:,i);
   else
    fi = 0;
 end
dw = dw - r(i)*fi;
end
      end
end
==============END wmerit.m=====================
==============dlag.m===========================
function dl = dlag(x,lambda)
[n,nx] = size(x);
[m,mr] = size(lambda);
[f,g] =fgv(x);
[vc,DC,m1,m] = constr(x);
dl = g - DC*(lambda);
end
=============END dlag.m========================
=============fgv6.m============================
function [f,g] = fgv(x)
[n, m]=size(x);
g=zeros(n,1);
f= (1.0 - x(1))^(2);
g(1)= 2.0*(x(1) - 1.0);
g(2)= 0.0;
end
============END fgv6.m========================
============constr6.m=========================
function [vc, DC, m1, m] = constr(x)
m = 1;
m1 = 1;
vc = zeros(m,1);
DC = zeros(2,m);
vc(1) = 10.0*(x(2) - x(1)^(2));
DC(1,1) = -20.0*x(1);
DC(2,1) = 10.0;
end
============END constr6.m===================
============fgv9.m==========================
function [f,g] = fgv(x)
[n, m]=size(x);
g=zeros(n,1);
f= sin(pi*x(1)/12.0)*cos(pi*x(2)/16.0);
g(1)= pi*cos(pi*x(1)/12.0)*cos(pi*x(2)/16.0)/12.0;
g(2)= -pi*sin(pi*x(1)/12.0)*sin(pi*x(2)/16.0)/16.0;
end
============END fgv9.m======================
============constr9.m=======================
function [vc, DC, m1, m] = constr(x)
m = 1;
m1 = 1;
vc = zeros(m,1);
DC = zeros(2,m);
vc(1) = 4.0*x(1) - 3.0*x(2);
DC(1,1) = 4.0;
DC(2,1) = -3.0;
end
===========END constr9.m==================
===========fgv14.m========================
function [f,g] = fgv(x)
[n, m]=size(x);
g=zeros(n,1);
f= (x(1) - 2)^(2) + (x(2) - 1)^(2);
g(1)= 2.0*(x(1) - 2);
g(2)= 2.0*(x(2) - 1);
end
==========END fgv14.m======================
==========constr14.m=======================
function [vc, DC, m1, m] = constr(x)
m = 2;
m1 = 1;
vc = zeros(m,1);
DC = zeros(2,2);
vc(1) = x(1) - 2.0*x(2) + 1;
vc(2) = -0.25*(x(1))^(2) - (x(2))^(2) + 1;
DC(1,1) = 1.0;
DC(2,1) = -2.0;
DC(1,2) = -0.5*x(1);
DC(2,2) = -2.0*x(2);
end
=========END constr14.m====================
=========fgv28.m===========================
function [f,g] = fgv(x)
[n, m]=size(x);
g=zeros(n,1);
f= (x(1) + x(2))^(2) + (x(2) + x(3))^(2);
g(1)= 2.0*(x(1) + x(2));
g(2)= 2.0*(x(1) + x(2)) + 2.0*(x(2) + x(3));
g(3)= 2.0*(x(2) + x(3));
end
===========END fgv28.m=====================
===========constr28.m======================
function [vc, DC, m1, m] = constr(x)
m = 1;
m1 = 1;
vc = zeros(m,1);
DC = zeros(3,1);
vc(1) = x(1) + 2.0*(x(2)) + 3.0*(x(3)) - 1;
DC(1,1) = 1.0;
DC(2,1) = 2.0;
DC(3,1) = 3.0;
end
===========END constr28.m====================
===========fgv32.m===========================
function [f,g] = fgv(x)
[n, m]=size(x);
g=zeros(n,1);
f= (x(1) + 3.0*x(2) + x(3))^(2) + 4.0*(x(1) - x(2))^(2);
g(1)= 2.0*(x(1) + 3.0*x(2) + x(3)) + 8.0*(x(1) - x(2));
g(2)= 6.0*(x(1) + 3.0*x(2) + x(3)) + 8.0*(x(2) - x(1));
g(3)= 2.0*(x(1) + 3.0*x(2) + x(3));
end
===========END fgv32.m======================
===========constr32.m=======================
function [vc, DC, m1, m] = constr(x)
m = 5;
m1 = 1;
vc = zeros(m,1);
DC = zeros(3,m);
vc(1) = -x(1) - (x(2)) - (x(3)) + 1;
vc(2) = 6.0*x(2) + 4.0*x(3) - (x(1))^(3) - 3.0;
vc(3) = x(1);
vc(4) = x(2);
vc(5) = x(3);
DC(1,1) = -1.0;
DC(2,1) = -1.0;
DC(3,1) = -1.0;
DC(1,2) = -3.0*(x(1))^(2);
DC(2,2) = 6.0;
DC(3,2) = 4.0;
DC(1,3) = 1.0;
DC(2,4) = 1.0;
DC(3,5) = 1.0;
end
=============END constr32.m=================
=============fgv39.m========================
function [f,g] = fgv(x)
[n, m]=size(x);
g=zeros(n,1);
f= -x(1);
g(1)= -1.0;
end
============END fgv39.m====================
============constr39.m=====================
function [vc, DC, m1, m] = constr(x)
m = 2;
m1 = 2;
vc = zeros(m,1);
DC = zeros(4,m);
vc(1) = -(x(1))^(3) + (x(2)) - (x(3))^(2);
vc(2) = (x(1))^(2) - x(2) - (x(4))^(2);
DC(1,1) = -3.0*(x(1))^(2);
DC(2,1) = 1.0;
DC(3,1) = -2.0*(x(3));
DC(1,2) = 2.0*x(1);
DC(2,2) = -1.0;
DC(4,2) = -2.0*x(4);
end
============END constr39.m================
