Netwon algorithm with backtracking is failing

7 次查看（过去 30 天）
Konstantinos
Konstantinos 2024-12-4
回答： Konstantinos 2024-12-6
Ran in:
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.
  2 个评论
Torsten
Torsten 2024-12-4
You forgot to define f (see above).
Please test your code for completeness before posting.
Konstantinos
Konstantinos 2024-12-4
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!

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-12-4
编辑：Torsten 2024-12-4
Ran in:
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
  2 个评论
Konstantinos
Konstantinos 2024-12-4
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
Torsten 2024-12-4
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.

请先登录，再进行评论。

更多回答（1 个）

Konstantinos
Konstantinos 2024-12-6
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

类别

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by