生成具有特征值的随机矩阵

Posted 2021-05-06

我正在做以下操作来生成具有特定范围内的特征值的随机矩阵：

function mat = randEig(dim, rReal)

    D=diff(rReal).*rand(dim,1)+rReal(1);
    P=rand(dim);
    mat=P*diag(D)/P;

end

但我也希望能够生成具有复杂（共轭）特征值的随机实矩阵。如何做到这一点？相似变换技巧将返回复杂矩阵。

编辑：好的，我设法通过搭载MATLAB的cdf2rdf函数（基本上是下面的第二个函数）来做到这一点。

function mat = randEig(dim, rangeEig, nComplex)

    if 2*nComplex > dim
         error('Cannot happen');
    end

    if nComplex
        cMat=diff(rangeEig).*rand(dim-2*nComplex,1)+rangeEig(1);
        for k=1:nComplex
            rpart=(diff(rangeEig).*rand(1,1)+rangeEig(1))*ones(2,1);
            ipart=(diff(rangeEig).*rand(1,1)+rangeEig(1))*i;
            ipart=[ipart; -ipart];
            cMat=[cMat; rpart+ipart];
        end
    else
        cMat=diff(rangeEig).*rand(dim,1)+rangeEig(1);
    end

    D=cMat;
    realDform = comp2rdf(diag(D));
    P=rand(dim);
    mat=P*realDform/P;
end


function dd = comp2rdf(d)
    i = find(imag(diag(d))');
    index = i(1:2:length(i));
    if isempty(index)
        dd=d;   
    else   
    if (max(index)==size(d,1)) | any(conj(d(index,index))~=d(index+1,index+1))
      error(message('Complex conjugacy not satisfied'));
    end
    j = sqrt(-1);
    t = eye(length(d));
    twobytwo = [1 1;j -j];
    for i=index
        t(i:i+1,i:i+1) = twobytwo;
    end 
       dd=t*d/t;
    end
end

但代码很丑陋，主要是多次调用rand的方式很烦人）。如果有人想发布一次调用rand的答案，并设法做到这一点，我一定会接受并赞成。

答案

我做了一个电话或两个电话：

function mat = randEig(dim, rangeEig, nComplex)

if 2*nComplex > dim
    error('Cannot happen');
end

if nComplex
    cMat=diff(rangeEig).*rand(2*nComplex,1)+rangeEig(1);
    cPart=cMat(1:nComplex)*i;
    cMat(1:nComplex)=[];
    cPart=upsample(cPart,2);
    cPart=cPart+circshift(-cPart,1);
    cMat=upsample(cMat,2);
    cMat=cMat+circshift(cMat,1);
    cMat=cMat+cPart;
    cMat=[diff(rangeEig).*rand(dim-2*nComplex,1)+rangeEig(1); cMat];
else
    cMat=diff(rangeEig).*rand(dim,1)+rangeEig(1);
end
D=cMat;
realDform = comp2rdf(diag(D));
P=rand(dim);
mat=P*realDform/P;
end


function dd = comp2rdf(d)
i = find(imag(diag(d))');
index = i(1:2:length(i));
if isempty(index)
    dd=d;
else
    if (max(index)==size(d,1)) | any(conj(d(index,index))~=d(index+1,index+1))
        error(message('Complex conjugacy not satisfied'));
    end
    j = sqrt(-1);
    t = eye(length(d));
    twobytwo = [1 1;j -j];
    for i=index
        t(i:i+1,i:i+1) = twobytwo;
    end
    dd=t*d/t;
end
end

如果有人可以拨打一个电话或更短/更优雅的代码，欢迎他们发布答案。