說明：

迭代法的收斂性和矩陣的條件數(shù)相關(guān)，條件數(shù)大于１K肯定不收斂，小于100肯定收斂

100--1000則要適當選擇截斷的小量，采用迭代法的另一種多參數(shù)調(diào)用方式

程序清單：

%%%%%%%%%%%%%%%%%%% 求最大特征值 %%%%%%%%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 估計矩陣的二范數(shù) %%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 結(jié)構(gòu)矩陣奇異值分解 %%%%%%%%%%%%%%%%% 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('錯誤，至少需要四個參數(shù)') 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個奇異值，因為矩陣階數(shù)不夠',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 矩陣乘以向量的運算，使用二倍擴充到循環(huán)矩陣用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函數(shù)，相當于一步迭代 % 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 %%

轉(zhuǎn)載于:https://www.cnblogs.com/xiao-cheng/archive/2011/12/09/2282473.html

總結(jié)

以上是生活随笔為你收集整理的Newton迭代法求解Toeplitz矩阵逆的程序的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。