Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss
Newton_radhson method
2 views
Skip to first unread message
zayed
Dec 31, 2011, 4:59:08 PM12/31/11
to
Hi,
I need to solve the following equation by using Newton-raphson method :
a. f(lambda)=-z'*inv(M)*inv(W)*Q*inv(M)*z,
W=inv(M)+lambda*Q
Q and W are complex square matrices ,also z and M are metrices .
I wrote a code to find the root by Newton's method ,but the value of the root is complex ,but it must be real.I am not sure a bout the derivative of f(x).
I have another form to the function f(x) ,but I don't know if it's suitable to be solved by Newton's method in matlab,the other form is:
m gamma_k*|x_k|^2
b. f(lambda)=sum ----------------------------------
k=1 (1+lambda(gamma_k)^2
where gamma_k is sub indices , gamma is eigenvalue
m is total number of eigenvalues of T=(M)^(1/3)*Q*(M)^(1/2)
x=U'*M^(-1/2)*z ,U is eigen vectors of T.
Is this equation suitable to be solved by Newton Raphson method.
The code for the form of (a) :
delta=1e-12;
epsilon=1e-12;
max1=500;
lambda=-1/(2*gamma);
for k=1:max1
% zeta=cos(theta);
% I=eye(n,n);
% Q=zeta*I-p*p';
W=inv(M)+lambda*Q;
y1=2*z'*inv(M)*inv(W)*Q.^2*inv(W)*inv(M)*z; % the derivative of y
y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
p1=lambda-y/y1;
err=abs(p1-lambda);
lambda=p1
if (err<delta)
break
end
k
err
end
I need to solve the following equation by using Newton-raphson method :
a. f(lambda)=-z'*inv(M)*inv(W)*Q*inv(M)*z,
W=inv(M)+lambda*Q
Q and W are complex square matrices ,also z and M are metrices .
I wrote a code to find the root by Newton's method ,but the value of the root is complex ,but it must be real.I am not sure a bout the derivative of f(x).
I have another form to the function f(x) ,but I don't know if it's suitable to be solved by Newton's method in matlab,the other form is:
m gamma_k*|x_k|^2
b. f(lambda)=sum ----------------------------------
k=1 (1+lambda(gamma_k)^2
where gamma_k is sub indices , gamma is eigenvalue
m is total number of eigenvalues of T=(M)^(1/3)*Q*(M)^(1/2)
x=U'*M^(-1/2)*z ,U is eigen vectors of T.
Is this equation suitable to be solved by Newton Raphson method.
The code for the form of (a) :
delta=1e-12;
epsilon=1e-12;
max1=500;
lambda=-1/(2*gamma);
for k=1:max1
% zeta=cos(theta);
% I=eye(n,n);
% Q=zeta*I-p*p';
W=inv(M)+lambda*Q;
y1=2*z'*inv(M)*inv(W)*Q.^2*inv(W)*inv(M)*z; % the derivative of y
y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
p1=lambda-y/y1;
err=abs(p1-lambda);
lambda=p1
if (err<delta)
break
end
k
err
end
Roger Stafford
Dec 31, 2011, 6:51:08 PM12/31/11
to
"zayed " <zaed...@Yahoo.com> wrote in message <jdo0jc$6g1$1...@newscl01ah.mathworks.com>...
- - - - - - - - - - -
Zayed, you say that "the root is complex, but it must be real". Why do you say that? In general the root of a complex equation would yield complex roots. That is in keeping with the philosophy that one should have two unknowns, namely the real and imaginary parts of the root, to solve two equations, namely the real and imaginary parts of the equation. Do you mean that in this particular case you expect the imaginary part of the root to be zero?
I don't quite follow your method of obtained the derivative of y with respect to lambda. For the sake of brevity I will use L to denote lambda. Also I will denote WI = inv(W) and MI = inv(M).
W = MI + L*Q, dW/dL = Q (Assuming M, MI, and Q are constant in lambda)
W*WI = I
dW/dL*WI + W*dWI/dL = dI/dL = 0
W*dWI/dL = -dW/dL*WI = -Q*WI
dWI/dL = -WI*Q*WI
y = -z'*MI*WI*Q*WI*MI*z
dy/dL = -z'*MI*(dWI/dL)*Q*WI*MI*z - z'*MI*WI*Q*(dWI/dL)*MI*z
= -z'*MI*(-WI*Q*WI)*Q*WI*MI*z - z'*MI*WI*Q*(-WI*Q*WI)*MI*z
= 2*z'*MI*WI*(Q*WI*Q)*WI*MI*z
whereas you obtained:
dy/dL = 2*z'*MI*WI*(Q.^2)*WI*MI*z
You have a Q.^2 in place of my Q*WI*Q . Are these supposed to be equal? It is not apparent to me. The first is presumably constant in L and the other variable. Could this be part of your trouble?
As to your question about form "b.", I see no reason why Newton-Raphson couldn't be used here. The derivative will be the sum of the derivatives of the m terms with respect to lambda and that looks easy to find if I understand your notation correctly.
Bear in mind that in using Newton-Raphson it is important to start your iteration with a good initial guess for the root. If there is more than one root, the root that is converged upon will depend on what initial value has been used. Also some badly selected initial values may cause the iteration to fail by, say, shooting off to infinity, or perhaps oscillating back and forth endlessly.
Roger Stafford
> I wrote a code to find the root by Newton's method ,but the value of the root is complex ,but it must be real.I am not sure a bout the derivative of f(x).
> .......
> W=inv(M)+lambda*Q;
> y1=2*z'*inv(M)*inv(W)*Q.^2*inv(W)*inv(M)*z; % the derivative of y
> y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
> .......
> y1=2*z'*inv(M)*inv(W)*Q.^2*inv(W)*inv(M)*z; % the derivative of y
> y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
- - - - - - - - - - -
Zayed, you say that "the root is complex, but it must be real". Why do you say that? In general the root of a complex equation would yield complex roots. That is in keeping with the philosophy that one should have two unknowns, namely the real and imaginary parts of the root, to solve two equations, namely the real and imaginary parts of the equation. Do you mean that in this particular case you expect the imaginary part of the root to be zero?
I don't quite follow your method of obtained the derivative of y with respect to lambda. For the sake of brevity I will use L to denote lambda. Also I will denote WI = inv(W) and MI = inv(M).
W = MI + L*Q, dW/dL = Q (Assuming M, MI, and Q are constant in lambda)
W*WI = I
dW/dL*WI + W*dWI/dL = dI/dL = 0
W*dWI/dL = -dW/dL*WI = -Q*WI
dWI/dL = -WI*Q*WI
y = -z'*MI*WI*Q*WI*MI*z
dy/dL = -z'*MI*(dWI/dL)*Q*WI*MI*z - z'*MI*WI*Q*(dWI/dL)*MI*z
= -z'*MI*(-WI*Q*WI)*Q*WI*MI*z - z'*MI*WI*Q*(-WI*Q*WI)*MI*z
= 2*z'*MI*WI*(Q*WI*Q)*WI*MI*z
whereas you obtained:
dy/dL = 2*z'*MI*WI*(Q.^2)*WI*MI*z
You have a Q.^2 in place of my Q*WI*Q . Are these supposed to be equal? It is not apparent to me. The first is presumably constant in L and the other variable. Could this be part of your trouble?
As to your question about form "b.", I see no reason why Newton-Raphson couldn't be used here. The derivative will be the sum of the derivatives of the m terms with respect to lambda and that looks easy to find if I understand your notation correctly.
Bear in mind that in using Newton-Raphson it is important to start your iteration with a good initial guess for the root. If there is more than one root, the root that is converged upon will depend on what initial value has been used. Also some badly selected initial values may cause the iteration to fail by, say, shooting off to infinity, or perhaps oscillating back and forth endlessly.
Roger Stafford
zayed
Jan 1, 2012, 12:44:08 PM1/1/12
to
"Roger Stafford" wrote in message <jdo75c$p09$1...@newscl01ah.mathworks.com>...
I would like to thank you Roger for your response . But I not sure a bout the initial value ,so I need help to plot the function in (a),so I can guess a good initial point.
Also I can't represent the form in (b) in matlab,so I can try to solve it.
Thanks
Also I can't represent the form in (b) in matlab,so I can try to solve it.
Thanks
Roger Stafford
Jan 1, 2012, 5:52:08 PM1/1/12
to
"zayed " <zaed...@Yahoo.com> wrote in message <jdq618$km9$1...@newscl01ah.mathworks.com>...
> > > ......
> > > y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
You have said that Q and W are complex-valued square matrices, and M would necessarily be square for you to speak of inv(M). Suppose you know three additional facts:
1. that M and Q are Hermitian (that is, M' = M and Q' = Q);
2. that lambda is assumed to be real;
3. and that z is a column vector (real or complex).
Then it immediately follows that inv(M), W, and inv(W) are also Hermitian. Also with z a column vector, then y must be a scalar. We can therefore write:
y' = (-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z)'
= -z'*inv(M)'*inv(W)'*Q'*inv(W)'*inv(M)'*z''
= -z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z
= y
which demonstrates that the scalar y must then be real-valued. (Note that I am using the prime symbol here, as in M', to denote matlab's conjugate transpose operator.)
Hence, given the above three conditions, you could conclude that y is always real. With complex values in Q (and possibly M) you may obtain very tiny imaginary parts for y due to round off error in actual computation, but you can then legitimately use only its real part by writing y = real(y). The question is, are these three above conditions really true for your problem? If so, that would permit you to make a simple plot of lambda versus y to determine approximately the real lambda value or values where y crosses over zero.
With respect to form (b) I think the same comments as above apply. If M and Q are Hermitian, then (using matlab's 'eig' function,) the eigenvalues of T will be real and so the indicated sum yields a real value for y. (Note: I assume you meant
T = (M)^(1/2)*Q*(M)^(1/2)
rather than M to the 1/3 power.) The 'eig' function also provides the eigenvectors of T which would enable you to find the 'x' values in this sum. I don't see what your difficulty is where you state "I can't represent the form in (b) in matlab." Why not?
Roger Stafford
> > "zayed " <zaed...@Yahoo.com> wrote in message <jdo0jc$6g1$1...@newscl01ah.mathworks.com>...
> > > W=inv(M)+lambda*Q;
> > > ......
> > > y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
> I would like to thank you Roger for your response . But I not sure a bout the initial value ,so I need help to plot the function in (a),so I can guess a good initial point.
> Also I can't represent the form in (b) in matlab,so I can try to solve it.
> Thanks
> Also I can't represent the form in (b) in matlab,so I can try to solve it.
> Thanks
- - - - - - - -
As to plotting your 'y' as a function of 'lambda' we very much need to know whether 'y' is strictly real-valued or not. That is because a "plot" of a complex-valued variable as a function of another complex-valued variable would require a four-dimensional display and matlab has no such direct display capability. Instead you would have to use a 3D surface plot for the separate real and imaginary parts of y as against the real and imaginary parts of lambda and determine where the two respective zero contour lines cross.
You have said that Q and W are complex-valued square matrices, and M would necessarily be square for you to speak of inv(M). Suppose you know three additional facts:
1. that M and Q are Hermitian (that is, M' = M and Q' = Q);
2. that lambda is assumed to be real;
3. and that z is a column vector (real or complex).
Then it immediately follows that inv(M), W, and inv(W) are also Hermitian. Also with z a column vector, then y must be a scalar. We can therefore write:
y' = (-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z)'
= -z'*inv(M)'*inv(W)'*Q'*inv(W)'*inv(M)'*z''
= -z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z
= y
which demonstrates that the scalar y must then be real-valued. (Note that I am using the prime symbol here, as in M', to denote matlab's conjugate transpose operator.)
Hence, given the above three conditions, you could conclude that y is always real. With complex values in Q (and possibly M) you may obtain very tiny imaginary parts for y due to round off error in actual computation, but you can then legitimately use only its real part by writing y = real(y). The question is, are these three above conditions really true for your problem? If so, that would permit you to make a simple plot of lambda versus y to determine approximately the real lambda value or values where y crosses over zero.
With respect to form (b) I think the same comments as above apply. If M and Q are Hermitian, then (using matlab's 'eig' function,) the eigenvalues of T will be real and so the indicated sum yields a real value for y. (Note: I assume you meant
T = (M)^(1/2)*Q*(M)^(1/2)
rather than M to the 1/3 power.) The 'eig' function also provides the eigenvectors of T which would enable you to find the 'x' values in this sum. I don't see what your difficulty is where you state "I can't represent the form in (b) in matlab." Why not?
Roger Stafford
zayed
Jan 1, 2012, 8:05:08 PM1/1/12
to
"Roger Stafford" wrote in message <jdqo2o$c07$1...@newscl01ah.mathworks.com>...
Thanks a lot Roger for your explanation it's useful.I derived the derivative of the function :
THE DERIVATIVE y1=2*z'*inv(M)*inv(W)*Q*inv(W)*Q*inv(W)*inv(M)*z;
THE FUNCTION y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
And the results is a complex numbers with a very small imaginary part ,also M,Q and W are Herimitian , is is column vector (real) ,Also the eigenvalues of T is complex with very little imaginary part.but the value of lambda is complex number with negative real part and very small imaginary part ,but lambda have to be real number and greater that zero ,Is this might be from the initial guess.
THE DERIVATIVE y1=2*z'*inv(M)*inv(W)*Q*inv(W)*Q*inv(W)*inv(M)*z;
THE FUNCTION y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
And the results is a complex numbers with a very small imaginary part ,also M,Q and W are Herimitian , is is column vector (real) ,Also the eigenvalues of T is complex with very little imaginary part.but the value of lambda is complex number with negative real part and very small imaginary part ,but lambda have to be real number and greater that zero ,Is this might be from the initial guess.
zayed
Jan 1, 2012, 8:07:08 PM1/1/12
to
Roger Stafford
Jan 1, 2012, 9:16:08 PM1/1/12
to
"zayed " <zaed...@Yahoo.com> wrote in message <jdqvs4$44b$1...@newscl01ah.mathworks.com>...
> Thanks a lot Roger for your explanation it's useful.I derived the derivative of the function :
> THE DERIVATIVE y1=2*z'*inv(M)*inv(W)*Q*inv(W)*Q*inv(W)*inv(M)*z;
> THE FUNCTION y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
>
> And the results is a complex numbers with a very small imaginary part ,also M,Q and W are Herimitian , is is column vector (real) ,Also the eigenvalues of T is complex with very little imaginary part.but the value of lambda is complex number with negative real part and very small imaginary part ,but lambda have to be real number and greater that zero ,Is this might be from the initial guess.
However, as I hinted to you in my last posting, you are necessarily doing your computation with less than perfect accuracy on a real-world computer. The binary floating point quantities you are using in these computations are only accurate out to about the sixteenth decimal place (52 bits) relative accuracy. We have no way of representing all possible numbers with perfect accuracy because that would require infinitely many bits.
Given this inaccuracy and the fact that complex quantities are involved, it is inevitable that you will obtain small imaginary parts to y, y1, and the eigenvalues of T. Their real parts will also be in error by similarly small amounts. This cannot be avoided. The sizes of these erroneous imaginary parts and the errors in the real parts will depend on the kinds of round off errors that are occurring, but with typical situations I would expect them to be roughly like a relative value of 1e-12 or thereabouts in relation to the sizes of the real parts, perhaps somewhat larger if your matrices are not well conditioned. You might try experimenting with the Symbolic Toolbox using extremely high accuracy, say 100 digits, to get a feel for the sizes of these errors.
If you are satisfied that these non-zero imaginary parts are just due to round off errors, you should simple discard them and use only the real parts. As I said, do y = real(y) and similarly with the eigenvalues of T.
Roger Stafford
> "Roger Stafford" wrote in message <jdqo2o$c07$1...@newscl01ah.mathworks.com>...
> > Hence, given the above three conditions, you could conclude that y is always real. With complex values in Q (and possibly M) you may obtain very tiny imaginary parts for y due to round off error in actual computation, but you can then legitimately use only its real part by writing y = real(y). .....
> Thanks a lot Roger for your explanation it's useful.I derived the derivative of the function :
> THE DERIVATIVE y1=2*z'*inv(M)*inv(W)*Q*inv(W)*Q*inv(W)*inv(M)*z;
> THE FUNCTION y=-z'*inv(M)*inv(W)*Q*inv(W)*inv(M)*z;
>
> And the results is a complex numbers with a very small imaginary part ,also M,Q and W are Herimitian , is is column vector (real) ,Also the eigenvalues of T is complex with very little imaginary part.but the value of lambda is complex number with negative real part and very small imaginary part ,but lambda have to be real number and greater that zero ,Is this might be from the initial guess.
- - - - - - - -
With perfect computation and perfectly Hermitian matrices along with real initial values of y, the computed values of y and its derivative, y1, as well as the eigenvalues of T would be real-valued. That is true in the ideal error-free mathematical world.
However, as I hinted to you in my last posting, you are necessarily doing your computation with less than perfect accuracy on a real-world computer. The binary floating point quantities you are using in these computations are only accurate out to about the sixteenth decimal place (52 bits) relative accuracy. We have no way of representing all possible numbers with perfect accuracy because that would require infinitely many bits.
Given this inaccuracy and the fact that complex quantities are involved, it is inevitable that you will obtain small imaginary parts to y, y1, and the eigenvalues of T. Their real parts will also be in error by similarly small amounts. This cannot be avoided. The sizes of these erroneous imaginary parts and the errors in the real parts will depend on the kinds of round off errors that are occurring, but with typical situations I would expect them to be roughly like a relative value of 1e-12 or thereabouts in relation to the sizes of the real parts, perhaps somewhat larger if your matrices are not well conditioned. You might try experimenting with the Symbolic Toolbox using extremely high accuracy, say 100 digits, to get a feel for the sizes of these errors.
If you are satisfied that these non-zero imaginary parts are just due to round off errors, you should simple discard them and use only the real parts. As I said, do y = real(y) and similarly with the eigenvalues of T.
Roger Stafford
zayed
Jan 1, 2012, 10:17:08 PM1/1/12
to
"Roger Stafford" wrote in message <jdr418$fbi$1...@newscl01ah.mathworks.com>...
The only problem is that the real part of lambda is negative ,and it should be positive real number .Is the reason round off errors.
Roger Stafford
Jan 1, 2012, 11:15:09 PM1/1/12
to
"zayed " <zaed...@Yahoo.com> wrote in message <jdr7jk$p8a$1...@newscl01ah.mathworks.com>...
I would recommend pursuing that plotting action you mentioned. That might tell you whether negative values of a real lambda can give rise to zero values of y.
I would also highly recommend following up on my suggestion to try directly computing y (not doing Newton-Raphson convergence) using the Symbolic Toolbox for the M and Q values you are using (make sure they are symbolically precisely Hermitian) and for various real values of lambda in the neighborhood of your supposed solutions. Matlab has a symbolic form of matrix inverse and matrix multiplication which would allow extremely high accuracy computation. To convert the results to numerical form you would need to use the 'vpa' function. You can set whatever accuracy you want for this. That should tell you whether or not your round off errors were responsible. However expect low speeds with this kind of computation.
Roger Stafford
> The only problem is that the real part of lambda is negative ,and it should be positive real number .Is the reason round off errors.
- - - - - - - -
It is impossible for me to answer that question without knowing more about your M and Q matrices and the z vector. It might be round off errors or it might possibly be multiple possible solutions for lambda with some of them being negative values. Have you checked that the solution you obtained satisfies y = 0 with sufficient accuracy?
I would recommend pursuing that plotting action you mentioned. That might tell you whether negative values of a real lambda can give rise to zero values of y.
I would also highly recommend following up on my suggestion to try directly computing y (not doing Newton-Raphson convergence) using the Symbolic Toolbox for the M and Q values you are using (make sure they are symbolically precisely Hermitian) and for various real values of lambda in the neighborhood of your supposed solutions. Matlab has a symbolic form of matrix inverse and matrix multiplication which would allow extremely high accuracy computation. To convert the results to numerical form you would need to use the 'vpa' function. You can set whatever accuracy you want for this. That should tell you whether or not your round off errors were responsible. However expect low speeds with this kind of computation.
Roger Stafford
zayed
Jan 2, 2012, 2:06:08 AM1/2/12
to
"Roger Stafford" wrote in message <jdrb0d$596$1...@newscl01ah.mathworks.com>...
Should I do the plot inside Newton iteration or in seperate action, because I tried to plot it in side newton iteration ,it's difficult to decide a bout the correct zero point .According to symbolic too box ,I am not familiar with this tool box ,so if it's possible I need more explanation in details
Roger Stafford
Jan 2, 2012, 2:47:08 AM1/2/12
to
"zayed " <zaed...@Yahoo.com> wrote in message <jdrl10$3fs$1...@newscl01ah.mathworks.com>...
For the purposes of making a plot, forget about the Newton-Raphson procedure for the moment. Just use a range of real lambda values, both negative and positive, which you think may cover all lambda values of possible interest and make a plot of y versus lambda. Or perhaps make several plots with gradually narrowing ranges and consequently higher resolution for areas where y is near zero. You won't need the derivative of y for this. For this purpose discard all the tiny imaginary parts of y. (However, first check that these imaginary values are indeed very small.) In this plot or these plots you will hopefully see the curve crossing zero at one or more points. This will give you an idea of the values of lambda that your Newton-Raphson procedure can later seek with much greater accuracy and it will also give you an educated guess about the appropriate initial starting point or
points to use there.
I cannot ease the path of understanding the Symbolic Toolbox. It will require some considerable work on your part reading the documentation. We all had to go through this learning process on that toolbox. However the effort will be well worth your time for a problem that may require great accuracy of computation to get results that you trust. In the past I have many times made good use of this facility. It is a part of matlab I value highly.
Roger Stafford
> Should I do the plot inside Newton iteration or in seperate action, because I tried to plot it in side newton iteration ,it's difficult to decide a bout the correct zero point .According to symbolic too box ,I am not familiar with this tool box ,so if it's possible I need more explanation in details
- - - - - - - - - - -
For the purposes of making a plot, forget about the Newton-Raphson procedure for the moment. Just use a range of real lambda values, both negative and positive, which you think may cover all lambda values of possible interest and make a plot of y versus lambda. Or perhaps make several plots with gradually narrowing ranges and consequently higher resolution for areas where y is near zero. You won't need the derivative of y for this. For this purpose discard all the tiny imaginary parts of y. (However, first check that these imaginary values are indeed very small.) In this plot or these plots you will hopefully see the curve crossing zero at one or more points. This will give you an idea of the values of lambda that your Newton-Raphson procedure can later seek with much greater accuracy and it will also give you an educated guess about the appropriate initial starting point or
points to use there.
I cannot ease the path of understanding the Symbolic Toolbox. It will require some considerable work on your part reading the documentation. We all had to go through this learning process on that toolbox. However the effort will be well worth your time for a problem that may require great accuracy of computation to get results that you trust. In the past I have many times made good use of this facility. It is a part of matlab I value highly.
Roger Stafford
0 new messages