Netwon algorithm with backtracking is failing

Hello everyone,

I am implementing the Newton method with backtracking line search, but I am encountering a problem with the dimensions of a matrix after a few iterations, and I can’t figure out why.

The algorithm I am asked to implement is as follows:

And the function that i am working with is the following:

.

Below is the code I am using:

n = 2; m = 20;               % Dimensions
c = randn(n, 1);               % Random vector c
A = randn(m, n);               % Random matrix A
b = abs(randn(m, 1))+1;          % Ensure b > 0

% Define the function
f = @(x)  (c' * x - sum(log(b - A * x)));

alpha = 0.3;
beta = 0.7;
iter_bt = 0;

%Newton method 

x_init = zeros(n,1);% Create a n by 1 dimensional vector

x_n = x_init;

thresh = 10^(-6);

f_gradient = @(x) c'- (A' * (1 ./ (b - A * x)));

F_newton = f(x_n); % Save function values for plotting
X_newton = x_n; % Save solutions for plotting

solved_newton = false;

while ~solved_newton


vec = 1 ./ ((b - A * x_n).^2);


diagonal = diag(vec);

f_hessian = @(x) A' * diagonal * A;

    % Compute gradient and Hessian
    grad = f_gradient(x_n);
    hess = f_hessian(x_n);

    % Compute Newton direction
    dxn_t = -hess \ grad;

    % Compute Newton decrement
    lambda_square = grad' * (hess \ grad); % Cancel the - sign in the dxn_t 
    
    % Check termination condition
    if (lambda_square / 2) <= thresh
        solved_newton = true;
        continue;
    end

    % Backtracking line search
    tau = backtracking_line_search(f, grad, x_n, dxn_t, alpha, beta, A, b);

    % Update iterate
    x_n = x_n + tau * dxn_t;
end

Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.

Error in solution>@(x)A'*diagonal*A (line 36)
f_hessian = @(x) A' * diagonal * A;

fprintf('Newton method converged to solution.\n');

function [tau] = backtracking_line_search(f, grad_f, x, direction, alpha, beta, A, b)
    % Backtracking line search
    % Inputs:
    %   f: function handle for objective function
    %   grad_f: gradient of f at current x
    %   x: current point
    %   direction: search direction (e.g., -grad for gradient descent)
    %   alpha: condition parameter (e.g., 0.3)
    %   beta: condition parameter (e.g., 0.7)
    %   A, b: constraints for feasibility check (b - A*x > 0)
    % Output:
    %   tau: step size satisfying the conditions
    %
    %Note:
    %    This implementation of the backtracking line search differs 
    %    from the original due to the necessity of the feasibility check. 
    %    Additionally, the order of the algorithm is slightly altered, 
    %    though this should not present an issue.
    
    tau = 1; % Initialize step size
    while true
        candidate_x = x + tau * direction; % Candidate point
        if all(b - A * candidate_x > 0) % Feasibility check
            if f(candidate_x) <= f(x) + alpha * tau * grad_f' * direction
                break; % Feasible and satisfies all the conditions
            end
        end
        tau = beta * tau; % Reduce step size
    end
end

After a few iterations of the while loop, MATLAB throws an error related to dimension mismatch during matrix operations. The issue seems to arise when constructing the Hessian matrix using f_hessian. I suspect something might be wrong with the way diag(vec) or A' * diag(vec) * A is computed.

Any insights, suggestions, or debugging tips would be greatly appreciated! Thank you in advance for your help.

[Torsten]

You forgot to define f (see above).

Please test your code for completeness before posting.

[Konstantinos]

I’ll update my question with the complete definition of f shortly to avoid any confusion. Thanks again for taking the time to review my post!

Answer 1

[Torsten]

rng('default')
n = 2; m = 20;
c = sym('c',[n 1]);
A = sym('A',[m n]);
b = sym('b',[m 1]);
x = sym('x',[n 1]);
f = c.'*x - sum(log(b-A*x));
grad = gradient(f,x)
hess = hessian(f,x);
f = matlabFunction(f,"Vars",{x,A,b,c});
f_gradient = matlabFunction(grad,"Vars",{x,A,b,c});
f_hessian = matlabFunction(hess,"Vars",{x,A,b,c});


A = randn(m, n);               % Random matrix A
b = abs(randn(m, 1))+1;          % Ensure b > 0
c = randn(n, 1);               % Random vector c

alpha = 0.3;
beta = 0.7;
iter_bt = 0;
%Newton method 
x_init = zeros(n,1); % Create a n by 1 dimensional vector
x_n = x_init

x_n = 2×1

0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

thresh = 10^(-6);

F_newton = f(x_n,A,b,c) % Save function values for plotting

F_newton = -8.6571

X_newton = x_n; % Save solutions for plotting
solved_newton = false;
while ~solved_newton

    % Compute gradient and Hessian
    grad = f_gradient(x_n,A,b,c);
    hess = f_hessian(x_n,A,b,c);
    % Compute Newton direction
    dxn_t = -hess \ grad;
    % Compute Newton decrement
    lambda_square = -grad' * dxn_t; % Cancel the - sign in the dxn_t 
    
    % Check termination condition
    if (lambda_square / 2) <= thresh
        solved_newton = true;
        continue;
    end
    % Backtracking line search
    tau = backtracking_line_search(f, grad, x_n, dxn_t, alpha, beta, A, b, c);
    % Update iterate
    x_n = x_n + tau * dxn_t
    f(x_n,A,b,c)
end

x_n = 2×1

-0.3456 -0.1127

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = -10.0602

x_n = 2×1

-0.4058 -0.1641

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = -10.0936

x_n = 2×1

-0.4020 -0.1613

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = -10.0938

fprintf('Newton method converged to solution.\n');

Newton method converged to solution.

function [tau] = backtracking_line_search(f, grad_f, x, direction, alpha, beta, A, b, c)
    % Backtracking line search
    % Inputs:
    %   f: function handle for objective function
    %   grad_f: gradient of f at current x
    %   x: current point
    %   direction: search direction (e.g., -grad for gradient descent)
    %   alpha: condition parameter (e.g., 0.3)
    %   beta: condition parameter (e.g., 0.7)
    %   A, b: constraints for feasibility check (b - A*x > 0)
    % Output:
    %   tau: step size satisfying the conditions
    %
    %Note:
    %    This implementation of the backtracking line search differs 
    %    from the original due to the necessity of the feasibility check. 
    %    Additionally, the order of the algorithm is slightly altered, 
    %    though this should not present an issue.
    
    tau = 1; % Initialize step size
    while true
        candidate_x = x + tau * direction; % Candidate point
        if all(b - A * candidate_x > 0) % Feasibility check
            if f(candidate_x,A,b,c) <= f(x,A,b,c) + alpha * tau * grad_f' * direction
                break; % Feasible and satisfies all the conditions
            end
        end
        tau = beta * tau; % Reduce step size
    end
end

[Konstantinos]

Is it necessary to use symbolic matrices as part of the solution?

In any case, thank you so much for your help! I will test the solution later and make sure to accept it. Thanks again for your support!

[Torsten]

Is it necessary to use symbolic matrices as part of the solution?

No. But it's the simplest way to get gradient and Hessian without errors.

If this was just a test problem and the dimensions of your "real" problem (m,n) are much larger, you should think about where you made a mistake in computing gradient and/or Hessian in order to avoid symbolic preprocessing.

Answer 2

[Konstantinos]

For completeness, here is a solution without the use of symbolic matrixes.

clearvars; clc; close all;

% Problem parameters
n = 2; m = 20;               % Dimensions
c = randn(n, 1);               % Random vector c
A = randn(m, n);               % Random matrix A
b = abs(randn(m, 1))+1;          % Ensure b > 0


alpha = 0.3;
beta = 0.7;
maxiters = 2000; %Use this in order to preallocate the arrays for faster execution time.
x_bt = zeros(n,1);
vals = zeros(1, maxiters);
thresh = 10^(-6);

iter_bt = 1; % Initialize iteration counter

while iter_bt <= maxiters %The algorithm will need less iterations than the max one.
    arg_log = b - A * x_bt;
    value = c' * x_bt - sum(log(arg_log));
    vals(iter_bt) = value;
    grad = c + A' * (1 ./ arg_log);
    hess = A'*diag(1./(arg_log.^2))*A;
    v = -hess\grad; %Prefer this way for inv(H).It provides more numerical stability
    lamba_square = grad' * (hess\grad);

    if (norm(grad) < thresh)
        break; 
    end
     

    %% Backtracking process
    
    fprime = grad' * v;
    tau = 1; % Initial step size
    darg_log = -(A * v) ./ arg_log;
    dc = c' * v;
    
    % Backtracking line search using while loop
    while true
        if any(b - A * (x_bt + tau * v) <= 0)
            tau = beta * tau;
        else
            newval = value + tau * dc - sum(log(1 + tau * darg_log));
            if newval > value + tau * alpha * fprime
                tau = beta * tau;
            else
                break;
            end
        end
    end
    
    % Update x
    x_bt = x_bt + tau * v;
    iter_bt = iter_bt + 1; % Increment iteration counter
end