0001 function eta = solve_along_line(M, point, x, y, g, Hy, sigma)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 im_tol = 1e-05;
0018
0019 inner = @(u, v) M.inner(point, u, v);
0020 rnorm = @(u) M.norm(point, u);
0021
0022 xx = inner(x, x);
0023 xy = inner(x, y);
0024 yy = inner(y, y);
0025 yHy = inner(y, Hy);
0026 const = inner(x, Hy) + inner(g, y);
0027
0028 func = @(a) a * const + 0.5 * a^2 * yHy + (sigma/3) * rnorm(M.lincomb(point, 1, x, a, y))^3;
0029
0030
0031
0032
0033
0034
0035 s2 = sigma * sigma;
0036 A = s2 * yy^3;
0037 B = 4 * s2 * xy * yy^2;
0038 C = 5 * s2 * xy^2 * yy + s2 * xx * yy^2 - yHy^2;
0039 D = 2 * s2 * xy * (xy^2 + xx * yy) - 2 * yHy * const;
0040 E = s2 * xx * xy^2 - const^2;
0041
0042 coeffs = [A, B, C, D, E];
0043 poly_roots = roots(coeffs);
0044 eta = 0;
0045 min_val = func(0);
0046 for root = poly_roots.'
0047 if root < 0 || abs(imag(root)) > im_tol
0048 continue;
0049 end
0050 rroot = real(root);
0051 root_val = func(rroot);
0052 if root_val < min_val
0053 eta = rroot;
0054 min_val = root_val;
0055 end
0056 end
0057
0058
0059
0060
0061
0062
0063 end