Quadratic Optimization using Newtons Method

I am writing Newtons method to solve an optimization problem. I have written the following code:

   function estimate = newtonMethod(f,df,x0,tolerance) %f & df are functions that contain OP & gradient respectively, x0 is initial guess, tolerance is 1e-12
   estimate = x0; % Initial guess 
    iterations = 40; % Maximum number of iterations
    for j = 1:iterations;
        dx = f(estimate) / df(estimate)';
        estimate = estimate - dx';
        % Stop criterion:
        if abs(dx) < tolerance
            return
        end
    end
end

However, the vector 'u' doesn't seem to be a good solution. My question is: Is this the correct way of solving this OP using Newtons Method? I would highly appreciate suggestions on this.

[Torsten]

Which OP ?

Best wishes

Torsten.

[rihab]

Its a Quadratic OP expressed as: f = sigma(k)norm(u' * A_k * u - b_k) ; with gradient df = sigma(k)4*(u'* A_k * u - b_k) * A_k * u (I have found the gradient by matrix calculus;A_k is symmetric,positive definite).

[Torsten]

Why a quadratic OP ? You take the norm of a quadratic in u - the result will have terms of order 4 in u, and a square root, too.

Best wishes

Torsten.

[rihab]

I see. Eventually I am interested in the square norm, and I am trying to solve it by Newtons Method (that approximates it to a quadratic problem). Somehow I see that the estimate generated by the code i have written (mentioned above) isn't a good one. Is this the correct way of implementing it in MATLAB?

[Torsten]

Why not use fminsearch again for control ?

Best wishes

Torsten.

[rihab]

I have already solved it using fminsearch which gave a very close estimate. I just want to see if Newtons method yields the same estimate or not so I wrote the code for it and wanted to know if its the correct way of implementing it or not :-)

[Torsten]

If you want to use Newton's method, you must apply it to the equation df/du = 0.

So you will have to calculate the Hessian of f: d^2f/du^2.

Best wishes

Torsten.

[rihab]

I see. So you mean instead of dx = f(estimate) / df(estimate); i should write dx = df(estimate) / dfHessian(estimate); (where df is the gradient and dfHessian is hessian of f)

[Torsten]

Yes.

[rihab]

thank you so much Torsten :-)