In [28]:

    

function newtmin_Toms566(p, maxIts, optTol)
    # obj is a function that evaluates the objective value, 
    #    gradient and hessian a x.
    # x0 is the starting point.
    # maxIts = maximum interations
    # optTol is the optimality tolerance on the gradient
    
    its = 1;
    Tol = 0;
    err = 0;
    fault_count = 0;
    
    # calculate the gradient from the objective function.  In our case, 
    # these values will be calls to the functions in Toms566
    x_k0 = p.x0;
    g_k0 = p.grd(x_k0);
    h_k0 = p.hes(x_k0);

    # Start the Loop 
    for i=1:maxIts
        # check for positive definiteness and full rank
        if( (det(h_k0)<0) && (rank(h_k0)==p.n) )
            # Spectral Decomposition of the Hessian
            H_diag,H_eigMat = eig(h_k0);
            # Replace all negative eigenvalues with their absolute value
                for j=1:length(H_diag)
                    if H_diag[j]<0
                        H_diag[j]=(abs(H_diag[j]));
                    end
                end
            h_k0 = H_eigMat*diagm(H_diag)*H_eigMat';
            fault_count+=1;

        end
        
        # Singular Hessian, exit the function
        if((det(h_k0)<0) && !(rank(h_k0)==p.n) )  
                err = 1;
            break;
        end
        
        #choose the Hessian/Modified Hessian
        B = h_k0;
        
        step_k0 = B\g_k0;
                
        #calculate the step length ----Armijo Backtracking 
        alpha = 1;
        while( p.obj(x_k0-alpha*step_k0) > ( (p.obj(x_k0)) ) )
            alpha=.5*alpha;
            if alpha < ( 1e-3)
                break;
            end
        end
        
        x_k1 = x_k0 - alpha*(step_k0);  
        x_k0 = x_k1;
        its+=1;
 
        g_k0 = p.grd(x_k0);
        h_k0 = p.hes(x_k0);
        
        Tol = norm(g_k0);
        if Tol<optTol
            break;
        end
    end
    x = x_k0;
    f = p.obj(x);
 
    return(f,its,Tol,h_k0,err,fault_count,x,alpha);
end


    

    
Out[28]:





newtmin_Toms566 (generic function with 1 method)

In [8]:

    

using Toms566


    

In [42]:

    

p=Problem(5)
n=p.n
f,its,Tol,h_k0,err,fault_count,x,alpha = newtmin_Toms566(p,30000,1e-10)


    

    
Out[42]:





(1.5759345474536323e-28,12,7.266682824072423e-15,
3x3 Array{Float64,2}:
  1.75192    -0.0298636    2.93685 
 -0.0298636   0.00475359  -0.113362
  2.93685    -0.113362     6.12801 ,

0,2,[0.999999999999949,10.000000000000414,1.0000000000000322],1)

In [44]:

    

A = convert(Array{Float64,2},readcsv("binary.csv")[2:end,:][:,1:3]);
using ForwardDiff

function L(x,u)
    a = x[1:2]; β = x[3];
    uᵢ = u[2:3]; yᵢ = u[1]
    return -yᵢ*(vecdot(a,uᵢ) + β) + log(1 + exp(vecdot(a,uᵢ) + β))
end;

G(x,ui) = ForwardDiff.gradient(x -> L(x,ui))(x);
H(x,ui) = ForwardDiff.hessian(x -> L(x,ui))(x);
m = size(A,1);
H(x) = sum([H(x, A[i,:]) for i in 1:m]);
L(x) = sum([L(x, A[i,:]) for i in 1:m]);
G(x) = sum([G(x, A[i,:]) for i in 1:m]);


    

In [26]:

    

function newtmin_LogReg(f, x_k0, maxIts, optTol)
    # obj is a function that evaluates the objective value, gradient and hessian at x.
    # x0 is the starting point.
    # maxIts = maximum interations
    # optTol is the optimality tolerance on the gradient
    
    its = 1;
    Tol = 0;
    err = 0;
    fault_count = 0;
       
    # intial values of the function
    x = x_k0;
    
    # calculate the gradient, hessian and objective at x 
    (f_k0, g_k0, h_k0) = f(x);
      
    # Start the Loop 
    for i=1:maxIts
        # check for positive definiteness and full rank
        if( (det(h_k0)<0) && (rank(h_k0)==3) )
            # Spectral Decomposition of the Hessian
            H_diag,H_eigMat = eig(h_k0);
            # Replace all negative eigenvalues with their absolute value
                for j=1:length(H_diag)
                    if H_diag[j]<0
                        H_diag[j]=(abs(H_diag[j]));
                    end
                end
            h_k0 = H_eigMat*diagm(H_diag)*H_eigMat';
            fault_count+=1;

       end
        # Singular Hessian, exit the function
        if((det(h_k0)<0) && !(rank(h_k0)==3) )  
                err = 1;
            break;
        end
        
        #choose the Hessian/Modified Hessian
        B = h_k0;
        
        step_k0 = B\g_k0;
                
        #calculate the step length ----Armijo Backtracking 
        alpha = 1;
        x=x_k0-alpha*step_k0;
        (f_k1, g_k1, h_k1) = f(x);
        while( f_k1 > f_k0 )
            alpha=.5*alpha;
            x=x_k0-alpha*step_k0;
            (f_k1, g_k1, h_k1) = f(x);
            if alpha < ( 1e-2)
                break;
            end
        end
        
        x_k1 = x_k0 - alpha*(step_k0);  
        x_k0 = x_k1;
        its+=1;
 
        x=x_k0;
        (f_k0, g_k0, h_k0) = f(x);
        
        Tol = norm(g_k0);
        if Tol<optTol
            break;
        end
    end
    x_out = x_k0;
    f_out = f_k0;
 
    return(f_out,its,Tol,h_k0,err,fault_count,x_out);
end


    

    
Out[26]:





newtmin_LogReg (generic function with 1 method)

In [43]:

    

f_out,its,Tol,h_k0,err,fault_count,x_out = newtmin_LogReg( x -> (L(x), G(x), H(x)) , zeros(3), 100, 1e-6)


    

    
Out[43]:





(240.17199084241437,6,3.242144089428076e-12,
3x3 Array{Float64,2}:
     3.11768e7    1.72965e5  49888.3   
     1.72965e5  989.84         284.17  
 49888.3        284.17          82.5031,

0,0,[0.0026906835959642976,0.7546868559629449,-4.949378062622567])

In [ ]: