function eta = solve_along_line(M, point, x, y, g, Hy, sigma)
% Minimize the function h(eta) = f(x + eta*y) where
%     f(s) = <s, H[s]> + <g, s> + sigma * ||s||^3.
%
% Inputs: A manifold M, a point on the manifold, vectors x, y, g, and H[y]
%         on the tangent space T_(point)M, and a constant sigma.
%
% Outputs: minimizer eta if eta is positive real; otherwise returns eta = 0

% This file is part of Manopt: www.manopt.org.
% Original authors: May 2, 2019,
%    Bryan Zhu, Nicolas Boumal.
% Contributors:
% Change log: 
    
    % Magnitude tolerance for imaginary part of roots.
    im_tol = 1e-05;
    
    inner = @(u, v) M.inner(point, u, v);
    rnorm = @(u) M.norm(point, u);

    xx = inner(x, x);
    xy = inner(x, y);
    yy = inner(y, y);
    yHy = inner(y, Hy);
    const = inner(x, Hy) + inner(g, y);
    
    func = @(a) a * const + 0.5 * a^2 * yHy + (sigma/3) * rnorm(M.lincomb(point, 1, x, a, y))^3;
    
    % upper_bound = Inf;
    % if bound_upper
    %     upper_bound = 1 / (ytHy/yty + sigma * rnorm(x));
    % end
    
    s2 = sigma * sigma;
    A = s2 * yy^3;
    B = 4 * s2 * xy * yy^2;
    C = 5 * s2 * xy^2 * yy + s2 * xx * yy^2 - yHy^2;
    D = 2 * s2 * xy * (xy^2 + xx * yy) - 2 * yHy * const;
    E = s2 * xx * xy^2 - const^2;
    
    coeffs = [A, B, C, D, E];
    poly_roots = roots(coeffs);    
    eta = 0;
    min_val = func(0);
    for root = poly_roots.'
        if root < 0 || abs(imag(root)) > im_tol
            continue;
        end
        rroot = real(root);
        root_val = func(rroot);
        if root_val < min_val
            eta = rroot;
            min_val = root_val;
        end
    end
    % if bound_upper 
    %     bound_val = func(upper_bound);
    %     if bound_val < min_val
    %         eta = upper_bound;
    %     end
    % end
end