说明：

迭代法的收敛性和矩阵的条件数相关，条件数大于１K肯定不收敛，小于100肯定收敛

100--1000则要适当选择截断的小量，采用迭代法的另一种多参数调用方式

程序清单：

%%%%%%%%%%%%%%%%%%% 求最大特征值 %%%%%%%%%%%%%%%%%%%%%%%%%%
CompEig.mfunction aerfa=CompEig(n)
%n=length(r);
q=eye(n,1);
tempLab=1;
lab=0;
max=100;
step=0;
while abs(tempLab-lab)>1e-10tempLab=lab;z=ToepMultA(q);z=ToepMultAt(z);z=ToepMultA(z);z=ToepMultAt(z);q=z/sqrt(z'*z);s=ToepMultA(q);s=ToepMultAt(s);s=ToepMultA(s);s=ToepMultAt(s);lab=q'*s;step=step+1;if step>=maxwarning('MATLAB:CompEig:notconverge', '%s%d%s', ...'Newton iteration did not converge for ', step, ' iterations with tolerance 1e-10');break;end
end
aerfa=lab;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  求矩阵的低秩分解形式 %%%%%%%%%%%%%%%%%%%%%%%
Deta.mfunction detaY=Deta(Y)
[m,n]=size(Y);
detaY1=zeros(m,n);
detaY2=detaY1;
detaY1(2:m,:)=Y(1:m-1,:);
detaY2(:,1:n-1)=Y(:,2:n);
detaY=detaY1-detaY2;%%%%%%%%%%%%%%%%%  矩阵乘向量   %%%%%%%%%%%%%%%%%%%%%%%%%function y=FunA(x)
y=ToepMultAt([0;x(1:end-1)]);
y=ToepMultA(y);
y=ToepMultAt(y);
x=ToepMultAt(x);
x=ToepMultA(x);
x=ToepMultAt(x);
y=[0;x(1:end-1)]-y;
endfunction y=FunAt(x)
y=[x(2:end);0];
y=ToepMultA(y);
y=ToepMultAt(y);
y=ToepMultA(y);
x=ToepMultA(x);
x=ToepMultAt(x);
x=ToepMultA(x);
y=y-[x(2:end);0];
endfunction y=FunY(c,r,x,aerfa)
m=length(c);
n=length(r);
cy0=c/aerfa;
ry0=r/aerfa;
top=-ry0(2:n)';
right=[ry0(n:-1:max(2,n-m+2));cy0(1:m-n)];
y=zeros(m,1);
y(1)=top*x(1:end-1);
y(2:m)=x(end)*right;
endfunction y=FunYt(c,r,x,aerfa)
m=length(c);
n=length(r);
cy0=c/aerfa;
ry0=r/aerfa;
top=-ry0(2:n)';
right=[ry0(n:-1:max(2,n-m+2));cy0(1:m-n)];
y=zeros(n,1);
y(end)=right'*x(2:end);
y(1:n-1)=x(1)*top;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  FFT 求解矩阵乘以向量  %%%%%%%%%%%%%%%%%%
function result=ToepMultA(y)
%
%   function for a Toeplitz Matrix to multiply a vector by FFT
%   A=toeplitz(c,r) is a toeplitz matrix%
global m;
global n;
global N;
global fftcc;
y=[y;zeros(N-n,1)];
result=ifft(fftcc.*fft(y));
result(m+1:N)=[];function result=ToepMultAt(y)
%
%   function for a Toeplitz Matrix to multiply a vector by FFT
%   A=toeplitz(c,r) is a toeplitz matrix
%
global m;
global n;
global N;
global fftrr;
y=[y;zeros(N-m,1)];
result=ifft(fftrr.*fft(y));
result(n+1:N)=[];function y=LarFunA(x)
% [m,n]=size(A)
% [n,m]=size(A')
% length(x)=m+n
% length(x1)=n
% length(x2)=m
%  [ 0  A' ]        x1          A'*x2
%           *           =
%  [ A  0  ]       x2          A*x1
global m;
global n;
y=[FunAt(x(m+1:m+n));FunA(x(1:m))];function y=LarFunY(x)
global m;
global n;
y=[FunYt(x(n+1:m+n));FunY(x(1:n))];
normFest
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    估计矩阵的二范数   %%%%%%%%%%%%%%%%%%%%
function [e,cnt] = normFest(Fun,FunT,m,tol)
%NORMEST Estimate the matrix 2-norm.
%   NORMEST(Fun,FunT,m) is an estimate of the 2-norm of the matrix A.
%   NORMEST(S,tol) uses relative error tol instead of 1.e-6.
%   [nrm,cnt] = NORMEST(..) also gives the number of iterations used.
%
%   This function is intended primarily for sparse matrices,
%   although it works correctly and may be useful for large, full
%   matrices as well.  Use NORMEST when your problem is large
%   enough that NORM takes too long to compute and an approximate
%   norm is acceptable.
%
%   INPUTS:
%   Fun, FunT:
%   the function to express matrix A as Fun=@(x) A*x and FunT=@(x) A'*x .
%   m: the number of rows of the matrix
%   tol: relative error
%   OUTPUTS:
%   e: the estimated norm of A
%   cnt: the number of iterations used%   Written by Xiao-Cheng, Dec,2011
%
%   Reference:
%   mfile of the function NORMESTif nargin < 4, tol = 1.e-6; end
maxiter = 100; % should never take this many iterations.
x = abs(FunT(ones(m,1)));
cnt = 0;
e = norm(x);
if e == 0, return, end
x = x/e;
e0 = 0;
while abs(e-e0) > tol*ee0 = e;Ax = Fun(x);if nnz(Ax) == 0Ax = rand(size(Ax));endx = FunT(Fun(x));normx = norm(x);e = normx/norm(Ax);x = x/normx;cnt = cnt+1;if cnt > maxiterwarning('MATLAB:normest:notconverge', '%s%d%s%g', ...'NORMFEST did not converge for ', maxiter, ' iterations with tolerance ', tol);break;end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  ODG 矩阵向量乘积 %%%%%%%%%%%%%%%%%%%%
function result = MultiplyByODG(c1,U,sigma,V,x,m,n)
%
%  compute Y*x by its odg expression
%  where trunc(Y)=L(c1)-\sum_{j=1}^{k}sigma(j)*L(u_j)*U(Z_n*v_j)
%
%   Input:
%        [c1,U,sigma,V] -- the odg of matrix Y
%         x -- as above
%        [m,n]  --  size of matrix Y
%   Output:
%        result -- the result Y*x%   Written by Xiao-Cheng, Nov, 2011
%
%   Reference:
%   Computing Moore-Penrose Inverses of Toeplitz matrices by Newton's Iteration
%   cc=[c1;zeros(m,1)];%%rr=[c1(1);zeros(m,1);c1(m:-1:2)];if m<=nxx=[x(1:m);zeros(m,1)];elsexx=[x(1:n);zeros(2*m-n,1)];end%length of result 2m ||---> mresult=ifft(fft(cc).*fft(xx));k=length(sigma);zz=zeros(2*m,1);xx=[x;zeros(n,1)];fftxx=fft(xx);for i=1:kv=V(:,i);u=U(:,i);%UU  2n-by-2n%LL  2m-by-2mccu=[zeros(n+1,1);v(n-1:-1:1)];%%rru=[0;v(1:n-1);zeros(n,1)];ccl=[u;zeros(m,1)];%%rrl=[u(1);zeros(m,1);u(m:-1:2)];yy=ifft(fft(ccu).*fftxx);        if m<n %L(u) m-by-m And U(Zn*v) m-by-nyy=[yy(1:m);zeros(m,1)];else%L(u) m-by-n And U(Zn*v) n-by-nyy=[yy(1:n);zeros(2*m-n,1)];         end%zz=zz+sigma(i)*ifft(fft(ccl).*fft(yy));zz=zz+sigma(i)*(fft(ccl).*fft(yy));endzz=ifft(zz);result=result(1:m)-zz(1:m);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = MultiplyTByODG(c1,U,sigma,V,x,m,n)
%
%  compute Y'*x by its odg expression
%  where trunc(Y)=L(c1)-\sum_{j=1}^{k}sigma(j)*L(u_j)*U(Z_n*v_j)
%        Y'=L'(c1)-\sum_{j=1}^{k}sigma(j)*U'(Z_n*v_j)*L'(u_j)
%
%   Input:
%        [c1,U,sigma,V] -- the odg of matrix Y
%         x -- as above length-m
%        [m,n]  --  size of matrix Y
%   Output:
%        result -- the result Y'*x%   Written by Xiao-Cheng, Nov, 2011
%
%   Reference:
%   Computing Moore-Penrose Inverses of Toeplitz matrices by Newton's Iteration
%   %   variables:
%   result: the result to return
%   cc: the expended column of the toeplitz matrix
%   xx: the expended vector of the vectorresult=zeros(n,1);cc=[c1(1);zeros(m,1);c1(m:-1:2)];xx=[x;zeros(m,1)];fftxx=fft(xx);tt=ifft(fft(cc).*fftxx);if m<nresult(1:m)=tt(1:m);elseresult=tt(1:n);endk=length(sigma);zz=zeros(2*n,1);for j=1:ku=U(:,j);v=V(:,j);ccl=[u(1);zeros(m,1);u(m:-1:2)];yy=ifft(fft(ccl).*fftxx);if m<nyy=[yy(1:m);zeros(2*n-m,1)];elseyy=[yy(1:n);zeros(n,1)];endccu=[0;v(1:n-1);zeros(n,1)];%zz=zz+sigma(j)*ifft(fft(ccu).*fft(yy));zz=zz+sigma(j)*(fft(ccu).*fft(yy));endzz=ifft(zz);result=result-zz(1:n);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  结构矩阵奇异值分解   %%%%%%%%%%%%%%%%%
function [U,S,V]=mysvds(varargin)
%MYSVDS this function is used to compute the svd of the large matrix where
%the matrix-vector multiply can be easily computed, such as the large
%sparse matrix or some structed matrices. The function only accepts the
%function expression of the matrix as parameters.
% svds for large matrix
% Input:
%       LFunA(varargin{1}): the function of the matrix, which is given
%               by a n-vector, returns another A*x
%       LFunAt(varargin{2}): the function of the transform matrix,which
%               is given a m-vector, returns another A'*x
%       [m,n](varargin{3,4}): the size of columns of the matrix A
%       k(varargin{5}): the number of eigenvalues to compute
% Output:
%       U,S,V: U*S*V'=\sum_{i=1}^{k}(\sigma_{i}u_{i}v_{i}^{A})
%   Written by Xiao-Cheng, 2010
%   Modified 23,Nov,2010
%   Modified  1,Dec,2011
%
%   Dependence:
%   normFest.m: the function to estimate the 2-norm of the matrix.
%
if nargin < 4error('错误，至少需要四个参数')
end
LFunA=varargin{1};
LFunAt=varargin{2};
m=varargin{3};
n=varargin{4};
q=max(m,n);
nA=normFest(LFunA,LFunAt,m);%the normest of A
M=min(m,n);
if nargin < 5% default to choose the first 6 eigenvaluessprintf('mysvds 将会给出前面6个最大的奇异值和奇异向量')k=min([M,6]);
elsek=min(M,varargin{5});if k==Msprintf('最多只能得到%d个奇异值，因为矩阵阶数不够',M)end
end
% compute the first 2k eigenvalues an eigenvectors of the
% matrix [zeros(n,n),T';T,zeros(m,m)]
LarFunA=@(x) [LFunAt(x(n+(1:m)));LFunA(x(1:n))];
opt.issym=1;
opt.isreal=1;
opt.disp=0;
[W,S]=eigs(LarFunA,m+n,2*k,'la',opt);
if ~isreal(S) || ~isreal(W)error('MYSVDS:ComplexValuesReturned',...['eigs([0 A''; A 0]) returned complex values -' ...' singular values of A cannot be computed in this way.'])
end
d=diag(S);[dum,ind]=sort(abs(d));
d=d(ind);
W=W(:,ind);dtol = q * nA * sqrt(eps);
uvtol = m * sqrt(sqrt(eps));VV = W(1:n,:)' * W(1:n,:);
dVV = diag(VV);
UU = W(n+(1:m),:)' * W(n+(1:m),:);
dUU = diag(UU);
indpos = find((d > dtol) & (abs(dUU-0.5) <= uvtol) & (abs(dVV-0.5) <= uvtol));
indpos = indpos(1:min(end,k));
V = sqrt(2) * W(1:n,indpos);
s = d(indpos);
U = sqrt(2) * W(n+(1:m),indpos);% sort the singular values in descending order
[s,ind] = sort(s);
s = s(end:-1:1);
U = U(:,ind(end:-1:1));
S = diag(s);
V = V(:,ind(end:-1:1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  一次ODG 过程 %%%%%%%%%%%%%%%%%%%%%%%%
function [yr,Ur,sigmar,Vr]=odg(y,Uy,sigmay,Vy,Ua,sigmaa,Va,epson)
%
% algorithm to compute the odg
%       A is a Toeplitz matrix A=toeplitz(c,r)
%   Input:
%        c -- the first column of the toeplitz matrix A
%        r -- the first row of the toeplitz matrix A,satisfy r(1)=c(1)
%        y -- the first column of the matrix Y
%        [Uy,sigmay,Vy] -- the SVD of the displacement operator of Y
%                          detaY=Uy*diag(sigmay)*Vy'=Z_m*Y-Y*Z_n
%        [Ua,sigmaa,Va] -- the SVD of the displacement operator of A'*A*A'
%                          detaA=Ua*diag(sigmaa)*Va'=Z_m*A'*A*A'-A'*A*A'*Z_n
%        epson -- for the eps-displacement rank
%   Output:
%       [yr,Ur,sigmar,Vr] -- the aodg of W or the odg of Y_{i+1}%   Written by Xiao-Cheng, Aug, 2010
%   Modefied by Xiao-Cheng, Nov, 2011%   Reference:
%   Computing Moore-Penrose Inverses of Toeplitz matrices by Newton's Iteration
%                                                          Yimin Wei
%
%   Dependence:
%           MutiplyAA.m -- Compute A'*A*A'*y
%           MutiplyY.m -- Compute Y*x
%           MutiplyYT.m -- Compute Y'*xglobal m;
global n;% Step 1
%y=2y-Y*A'*A*A'*y
x=ToepMultAt(y);
x=ToepMultA(x);
x=ToepMultAt(x);
yr=2*y-MultiplyByODG(y,Uy,sigmay,Vy,x,m,n);%Step2
% Y*A'*A*A'*Uy
ky=length(sigmay);
ka=length(sigmaa);
TempY=Uy;
for i=1:kyx=ToepMultAt(TempY(:,i));x=ToepMultA(x);x=ToepMultAt(x);TempY(:,i)=MultiplyByODG(y,Uy,sigmay,Vy,x,m,n);
end% Y*Ua
s=zeros(m,ka);
for i=1:kas(:,i)=MultiplyByODG(y,Uy,sigmay,Vy,Ua(:,i),m,n);
endUw=[Uy,TempY,s];
%Uw=[Uy,Y*A'*A*A'*Uy,Y*Ua];% Y'*Va
s=zeros(n,ka);
for i=1:kas(:,i)=MultiplyTByODG(y,Uy,sigmay,Vy,Va(:,i),m,n);
end% Y'*A*A'*A*Vy
t=zeros(n,ky);
for i=1:kyx=ToepMultA(Vy(:,i));x=ToepMultAt(x);x=ToepMultA(x);t(:,i)=MultiplyTByODG(y,Uy,sigmay,Vy,x,m,n);
endVw=[Vy,s,t];
%Vw=[Vy,Y'*Va,Y'*A*A'*A*Vy];
%\Sigma_w
% Sw=[[2*diag(sigmay),zeros(ky,ka),-diag(sigmay)]; ...
%     -[diag(sigmay),zeros(ky,ka+ky)];...
%     [zeros(ka,ky),-diag(sigmaa),zeros(ka,ky)]]; % finished but not needed
%Uw*Sw*Vw';%Step 3
[Q1,R1]=qr(Uw,0);
[Q2,R2]=qr(Vw,0);%Step 4
% R1=[Rk1,Rk2,Rh]
% R1*Sw=[(2Rk1-Rk2)S,-RhSa,Rk1S]...column weighted%%%% compute R1=R1*Sw
R11=2*R1(:,1:ky)-R1(:,ky+1:2*ky);
for j=1:kaR1(:,ky+j)=-sigmaa(j)*R1(:,2*ky+j);
end
for i=1:kyR1(:,ky+ka+i)=-sigmay(i)*R1(:,i);R1(:,i)=R11(:,i)*sigmay(i);
end
%%%%
%%% compute R3=R1*Sw*R2'
R3=sparse(R1)*sparse(R2');[U,S,V]=svds(R3,ka+2*ky);sigmar=diag(S);
epson=sigmar(1)*epson;
tempr=find([sigmar;0]<=epson);
r=tempr(1)-1;%Step 5
Ur=Q1*U(:,1:r);
Vr=Q2*V(:,1:r);
sigmar=sigmar(1:r);% % test
% % DUw=[Uy,Y*A'*A*A'*Uy,Y*Ua];
% % DVw=[Vy,Y'*Va,Y'*A*A'*A*Vy];
%%DSw=[2*diag(sigmay),zeros(ky,ka),-diag(sigmay);-diag(sigmay),zeros(ky,ka),zeros(ky,ky);zeros(ka,ky),-diag(sigmaa),zeros(ka,ky)];
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  Newton 迭代 %%%%%%%%%%%%%%%%%%%%%%%%%
% Modified Newton's Iteration for computing the P-inverse of the
% toeplitz matrix
% Input:
%       c=varargin{1}: the first column of the toeplitz matrix
%       r=varargin{2}: the first row of the toeplitz matrix
%       epsoni=varargin{3}: a parameter for the iteration algorithm
% Output:
%       [y,Uy,sigmay,Vy]: the ODG of the matrix pinv(A)
%       remark: we can compute the P-inverse of A from the ODG
%
% Dependence:
%       mysvds.m -- Compute the svd of matrix
%       FunA.m -- support a function handle to compute A*x
%       FunAt -- support a function handle to compute A'*x
%       CompEig.m -- compute the largest eigenvalue of a Toeplitz matrix
%       FunY.m -- support a function handle to compute Y0*x
%       FunYt.m -- support a function handle to compute Y0'*x
function [y,Uy,sigmay,Vy]=ModNewIter(varargin)
c=varargin{1};
r=varargin{2};
m=length(c);
n=length(r);
if nargin==3epsoni=varargin{3};
end
% Step 1
% Compute the odg of A'*A*A'
% A=toeplitz(c,r);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  AA=zeros(n,m);
%  AAA=zeros(m,m);
%  for i=1:m
%      y=ToepMultA(A(i,:)');
%      AA(:,i)=ToepMultAt(y);
%      AAA(:,i)=ToepMultA(AA(:,i));
%  end
%  DetaA=Deta(AA);
% [Ua,Sa,Va]=svds(DetaA,6);
% norm(Ua*Sa*Va'-DetaA)[Ua,Sa,Va]=mysvds(@FunA,@FunAt,n,m,6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigmaa=diag(Sa);
% Setp 2
% choose aerfa, Compute the odg of Y_0
% aerfa
aerfa=abs(CompEig(n));
cy0=c/aerfa;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ry0=r/aerfa;
% Y0=toeplitz(cy0,ry0);
% [Uy,Sy,Vy]=svds(Deta(Y0),2);
LFunY=@(x) FunY(c,r,x,aerfa);
LFunYt=@(x) FunYt(c,r,x,aerfa);
[Uy,Sy,Vy]=mysvds(LFunY,LFunYt,m,n,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%test
% theta=norm(A*pinv(A)-AAA/aerfa)
% condA=cond(A)
%test
sigmay=diag(Sy);
y=cy0;% Step 3
% determine epsoni, compute the odg of Y_{i+1} by odg()
Ri = 1;
step = 0;
tol=1e-8;
normA=normFest(@ToepMultA,@ToepMultAt,m);
while Ri/normA > tolz=MultiplyTByODG(y,Uy,sigmay,Vy,c,m,n);z=ToepMultA(z);z=ToepMultAt(z);x=ToepMultAt(c);x=MultiplyByODG(y,Uy,sigmay,Vy,x,m,n);x=ToepMultAt(x);e4=norm(x-z);x=ToepMultA(x);e1=norm(c-x);x=MultiplyByODG(y,Uy,sigmay,Vy,r,m,n);x=ToepMultAt(x);z=x;z=ToepMultA(z);z=ToepMultAt(z);z=MultiplyByODG(y,Uy,sigmay,Vy,z,m,n);z=ToepMultAt(z);e2=norm(x-z);x=ToepMultA(x);z=ToepMultA(r);z=MultiplyTByODG(y,Uy,sigmay,Vy,z,m,n);z=ToepMultA(z);e3=norm(z-x);Ri=max([e1,e2,e3,e4])if Ri>1e5error('The Newton iteration cannot be convergence, choose another eps')endif step>=40warning('MATLAB:ModNewIter:notconverge', '%s%d%s%g', ...'Newton iteration did not converge for ', step, ' iterations with tolerance ', tol);break;endstep=step+1;if nargin<3epsoni=Ri/normA^4;end[y,Uy,sigmay,Vy]=odg(y,Uy,sigmay,Vy,Ua,sigmaa,Va,epsoni);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  测试  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
clear global;
global m;
global n;
global fftcc;
global fftrr;
global N;
global r;
global c;%%
% q=4;
% n=256;
% m=n-2*q;
% N=2*max(m,n);
% c=zeros(m,1);
% r=zeros(n,1);
% for i=1:2*q+1
%     r(i)=1/(2*q+1)^2;
% end
% c(1)=1/(2*q+1)^2;
%%
n=8192;
m=8192;
N=2*max(m,n);
c=zeros(m,1);
r=zeros(n,1);
for i=1:m-1c(i)= 1/i;
end
c(1)=1;
c(m)=1;
% r=flipud(c);
for i=1:nr(i)=1/(n-i+1);
end
r(1)=1;
%%
% m=128;
% n=128;
% N=2*max(m,n);
% c=zeros(m,1);
% r=zeros(n,1);
% for i=1:m-1
%     c(i)= rand(1);
% end
% c(1)=1;
% c(m)=1;
% % r=flipud(c);
% for i=1:n
%     r(i)=rand(1);
% end
% r(1)=1;%%
rr=[r;zeros(N-(m+n-1),1);c(m:-1:2)];
cc=[c;zeros(N-(m+n-1),1);r(n:-1:2)];
fftcc=fft(cc);
fftrr=fft(rr);%A=toeplitz(c,r);
%normA=norm(A);format long;%%%%%%%%%%%%%    test function for odg     %%%%%%%%%%%%%%%%
% A=toeplitz(c,r);
% normA=normest(A);
% condA=cond(A);
%%
%测试toeplitz 矩阵乘以向量的运算，使用二倍扩充到循环矩阵用FFT计算
% x1=rand(m,1);
% x2=rand(m,1);
% y1=rand(n,1);
% y2=rand(n,1);
% ErrorAy1=ToepMultA(y1);
% ErrorAy2=ToepMultA(y2);
% ErrorAtx1=ToepMultAt(x1);
% ErrorAtx2=ToepMultAt(x2);
% disp('Error:A*y');
% disp(norm(ToepMultA(y1-y2)-(ErrorAy1-ErrorAy2))/norm(y1-y2));
% disp('Error:A''*y');
% disp(norm(ToepMultAt(x1-x2)-(ErrorAtx1-ErrorAtx2))/norm(x1-x2));%%
%测试矩阵乘以向量，使用odg表达式形式
% x=rand(n,1);
% DetaA=zeros(m,n);
% DetaA(1,1:n-1)=-r(2:n);
% if m<=n
%     DetaA(2:m,n)=r(n:-1:n-m+2);
% else
%     DetaA(2:m,n)=[r(n:-1:2);c(1:m-n)];
% end
% [U,S,V]=svds(DetaA, 2);
% sigma=diag(S);
% tic
% resultByOdg=MultiplyByODG(c,U,sigma,V,x,m,n);
% toc
% tic
% resulty=MutiplyY(c,U,sigma,V,x);
% toc
% result=A*x;
% disp('Error:odg multiply');
% disp(norm(resultByOdg-result));
% disp(norm(resultByOdg-resulty));
%
% x=rand(m,1);
% tic
% resultTByOdgT=MultiplyTByODG(c,U,sigma,V,x,m,n);
% toc
% tic
% resulty=MutiplyYT(c,U,sigma,V,x);
% toc
% DetaAT=zeros(n,m);
% DetaAT(1,1:m-1)=-c(2:m);
% if n<=m
%     DetaAT(2:n,m)=c(m:-1:m-n+2);
% else
%     DetaAT(2:n,m)=[c(m:-1:2);r(1:n-m)];
% end
% [U,S,V]=svds(DetaAT, 2);
% sigma=diag(S);
% resultTByOdg=MultiplyByODG(r,U,sigma,V,x,n,m);
%  resultT=A'*x;
% disp(norm(resultTByOdg-resultT));
% disp(norm(resultTByOdgT-resultT));
% disp(norm(resultTByOdgT-resulty));
%%%%%%%%%%%%%    test function for odg    %%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%  test odg   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 测试odg函数，相当于一步迭代
% A=toeplitz(c,r);
% Y=0.1*A;
% W=2*Y-Y*A'*A*A'*Y;
%
% DetaA=Deta(A'*A*A');
% DetaY=Deta(Y);
% DetaW=Deta(W);
%
% [Ua,Sa,Va]=svds(DetaA, 6);
% [U,S,V]=svds(DetaY,2);
% %odg function
% [yr,Ur,sigmar,Vr]=odg(Y(:,1),U,diag(S),V,Ua,diag(Sa),Va,1e-8);
%
% errorDetaW=norm(DetaW-Ur*diag(sigmar)*Vr');
% q=length(sigmar);
% [Uw,Sw,Vw]=svds(DetaW,q);
% w=W(:,1);
% errorDetaWFirstCol=norm(yr-w);
% errorDiagDetaW=norm(sigmar-diag(Sw));
%
% disp('Error ODG:');
% disp('DetaW Error');
% disp(errorDetaW);
% disp('W''s first column error');
% disp(errorDetaWFirstCol);
% disp('sigle value error');
% disp(errorDiagDetaW);
%%%%%%%%%%%%%%%   test odg     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%test msvds
% F=@(x) A*x;
% FT=@(x) A'*x;
% %normA=normest(A);
% tic
% [Um,Sm,Vm]=mysvds(@ToepMultA,@ToepMultAt,m,n);
% toc
% tic
% [Uc,Sc,Vc]=svds(A);
% toc
% norm(Um'*Um)
% norm(Um*Sm*Vm'-A)
% norm(Sm-Sc)
% tic
% normA1=normest(A);
% toc
% tic
% normA2=normFest(@ToepMultA,@ToepMultAt,m);
% toc
% norm(normA1-normA2)%%
% %%%%%%%%%%%%%%%   test Newton Iteration  %%%%%%%%%%%%%%
epson=1e-9;
% if condA>500
%     disp('The condition nember of A is too large, Newton iteration cannot be convergence');
% else
tic
[y,Uy,sigmay,Vy]=ModNewIter(c,r);
toc
% t=zeros(m,m);
% for i=1:m
%     t(:,i)=MultiplyByODG(y,Uy,sigmay,Vy,A(i,:)',m,n);
% end
% s=zeros(n,m);
% for i=1:m
%     s(:,i)=ToepMultAt(t(:,i));
% end
% disp('The condition number of A is');
% disp(condA);
% PInverseError=norm(pinv(A)-s)/normest(s);
% disp('The relative M-P inverse error measured by 2-norm is');
% disp(PInverseError);
% end
%%