Newton迭代法求解Toeplitz矩阵逆的程序
说明：
迭代法的收敛性和矩阵的条件数相关，条件数大于１K肯定不收敛，小于100肯定收敛
100--1000则要适当选择截断的小量，采用迭代法的另一种多参数调用方式
程序清单：
%%%%%%%%%%%%%%%%%%% 求最大特征值 %%%%%%%%%%%%%%%%%%%%%%%%%%
CompEig.m
function aerfa=CompEig(n)
%n=length(r);
q=eye(n,1);
tempLab=1;
lab=0;
max=100;
step=0;
while abs(tempLab-lab)>1e-10
tempLab=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>=max
warning('MATLAB:CompEig:notconverge', '%s%d%s', ...
'Newton iteration did not converge for ', step, ' iterations with tolerance 1e-10');
break;
end
end
aerfa=lab;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 求矩阵的低秩分解形式 %%%%%%%%%%%%%%%%%%%%%%%
Deta.m
function 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;
end
function 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];
end
function 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;
end
function 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 NORMEST
if 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*e
e0 = e;
Ax = Fun(x);
if nnz(Ax) == 0
Ax = rand(size(Ax));
end
x = FunT(Fun(x));
normx = norm(x);
e = normx/norm(Ax);
x = x/normx;
cnt = cnt+1;
if cnt > maxiter
warning('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<=n
xx=[x(1:m);zeros(m,1)];
else
xx=[x(1:n);zeros(2*m-n,1)];
end
%length of result 2m ||---> m
result=ifft(fft(cc).*fft(xx));
k=length(sigma);
zz=zeros(2*m,1);
xx=[x;zeros(n,1)];
fftxx=fft(xx);
for i=1:k
v=V(:,i);
u=U(:,i);
%UU 2n-by-2n
%LL 2m-by-2m
ccu=[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-n
yy=[yy(1:m);zeros(m,1)];
else
%L(u) m-by-n And U(Zn*v) n-by-n
yy=[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));
end
zz=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 vector
result=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<n
result(1:m)=tt(1:m);
else
result=tt(1:n);
end
k=length(sigma);
zz=zeros(2*n,1);
for j=1:k
u=U(:,j);
v=V(:,j);
ccl=[u(1);zeros(m,1);u(m:-1:2)];
yy=ifft(fft(ccl).*fftxx);
if m<n
yy=[yy(1:m);zeros(2*n-m,1)];
else
yy=[yy(1:n);zeros(n,1)];
end
ccu=[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));
end
zz=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 < 4
error('错误，至少需要四个参数')
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 eigenvalues
sprintf('mysvds 将会给出前面6个最大的奇异值和奇异向量')
k=min([M,6]);
else
k=min(M,varargin{5});
if k==M
sprintf('最多只能得到%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'*x
global 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:ky
x=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:ka
s(:,i)=MultiplyByODG(y,Uy,sigmay,Vy,Ua(:,i),m,n);
end
Uw=[Uy,TempY,s];
%Uw=[Uy,Y*A'*A*A'*Uy,Y*Ua];
% Y'*Va
s=zeros(n,ka);
for i=1:ka
s(:,i)=MultiplyTByODG(y,Uy,sigmay,Vy,Va(:,i),m,n);
end
% Y'*A*A'*A*Vy
t=zeros(n,ky);
for i=1:ky
x=ToepMultA(Vy(:,i));
x=ToepMultAt(x);
x=ToepMultA(x);
t(:,i)=MultiplyTByODG(y,Uy,sigmay,Vy,x,m,n);
end
Vw=[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:ka
R1(:,ky+j)=-sigmaa(j)*R1(:,2*ky+j);
end
for i=1:ky
R1(:,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==3
epsoni=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 > tol
z=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>1e5
error('The Newton iteration cannot be convergence, choose another eps')
end
if step>=40
warning('MATLAB:ModNewIter:notconverge', '%s%d%s%g', ...
'Newton iteration did not converge for ', step, ' iterations with tolerance ', tol);
break;
end
step=step+1;
if nargin<3
epsoni=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-1
c(i)= 1/i;
end
c(1)=1;
c(m)=1;
% r=flipud(c);
for i=1:n
r(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
%%
浙公网安备 33010602011771号