MATLAB Newton Raphson Method with a vector function

Hi! I am trying to use Newton Raphson Method to solve the following problem: Given a complex matrix $K$, let our vector variable be , where and . Also, my function is . I want to find a vector $p$ such that . To do that I use the following Newton Raphson algorithm , where the is a Jacobian (if Jacobian is singular, I think we can use a Moore-Penrose). I tried to modify the code here: https://www.mathworks.com/matlabcentral/answers/753414-matlab-newton-raphson-method-with-a-function-with-array-matrix-variables?s_tid=mwa_osa_a

Unfortunately, it is not working properly, even when . Can you please help me to fix it?

K=2.*rand(2,2)-1*ones(2,2)+i*(2.*rand(2,2) -1*ones(2,2));
K1=K+K';
K2=K'*K;
I=eye(2,2);

f = @(x1,x2,x3,x4) [K1*[x1 x2]'-x3*K2*[x1 x2]'-x4*[x1 x2]'; [x1 x2]*K2*[x1 x2]'-1; [x1 x2]*[x1 x2]'-1];
df = @(x1,x2,x3,x4) [K1-x3*K2-x4*I -K2*[x1 x2]' -[x1 x2]'; [x1 x2]*K2 0 0; [x1 x2] 0 0];

n=10
for i= 1:n
    x10 = 0.001;
    x20 = 0.001;
    x30 = 0.001;
    x40 = 0.001;
    y = newtonraphson(f, df, x10, x20, x30, x40, n);
    fprintf('%.4f  \n', y);
end
function y = newtonraphson(f, df, x10, x20, x30, x40, n)
rounding = 5*10^-(4);

for i=1:100
    f0 = f(x10,x20,x30,x40); 
    f0_der=df(x10,x20,x30,x40); 
    y=[x10 x20 x30 x40]'-f0_der\f0; 
    err=abs(y-[x10 x20 x30 x40]');
    
    if err<ones(4,1)*rounding 
        break;
    end
    x10=y(1);
    x20=y(2);
    x30=y(3);
    x40=y(4);
end

fprintf('The Root is : %f \n',y);
fprintf('No. of Iterations : %d\n',i);
end

Answer 1

[Torsten]

You have 7 equations for 4 unknowns. This won't give an exact solution, but only a solution in the least-squares sense.

K=2.*rand(2,2)-1*ones(2,2)+i*(2.*rand(2,2) -1*ones(2,2));
K1=K+K';
K2=K'*K;
x = lsqnonlin(@(x) fun(x(1),x(2),x(3),x(4),K1,K2),ones(4,1))

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

x = 4×1

0.0091 0.8964 1.3065 -0.0835

res = norm(fun(x(1),x(2),x(3),x(4),K1,K2))

res = 0.5506

function res = fun(x1,x2,x3,x4,K1,K2);
  res = zeros(7,1);
  res(1:2) = real(K1*[x1 x2]'-x3*K2*[x1 x2]'-x4*[x1 x2]');
  res(3:4) = imag(K1*[x1 x2]'-x3*K2*[x1 x2]'-x4*[x1 x2]');
  res(5) = real([x1 x2]*K2*[x1 x2]'-1);
  res(6) = imag([x1 x2]*K2*[x1 x2]'-1);
  res(7) = [x1 x2]*[x1 x2]'-1;
end

[Torsten]

That's an optimization problem, not a root-finding problem.

Use MATLAB's "fmincon" to solve.