对角化稀疏酉矩阵

Question

我必须收集稀疏酉矩阵的特征值。 基本上只有一个元素不同于零 行和列（这是一些马尔可夫过程的传递矩阵）。

我的问题是如何进行，最好的选择是什么 在所有功能套件中。我已经看到 eigs 可以提供帮助， 但我也看到必须选择初始向量。

Answer 1

以下代码最终定义了pdeig，它返回一个矩阵的特征值，该矩阵是一个pdmatrix，即一个排列和对角矩阵的乘积，或者换句话说，一个像问题描述的矩阵。也可以快速计算特征向量（它们有一个明确的公式）：

issquare(m) = all(x->x==size(m,1),size(m))
isunique(v) = v == unique(v)
permmatrix(sigma) = 
  [i==sigma[j] ? 1.0 : 0.0 for i=1:length(sigma),j=1:length(sigma)]
mat2perm(m) = [findfirst(m[:,i]) for i=1:size(m,1)]

function ispdmatrix(m)      # used to verify input matrix form
  (r,c,v) = findnz(m)
  return issquare(m) && isunique(r) && isunique(c)
end

function pdfact(m::Matrix)  # factor into permutation/dilation
  ispdmatrix(m) || error("input matrix must be a PD matrix")
  n = size(m,1)
  p = mat2perm(m)
  d = [p[i]>0 ? m[p[i],i] : zero(eltype(m)) for i=1:n]
  return (p,d)
end

# return eigenvalues from factored pdmatrix
function pdeig(p::Vector{Int},d::Vector)
  n = length(p)
  active = trues(n)
  eigv = Vector{Complex{eltype(d)}}(0)
  for i=1:n
    if !active[i]
      continue
    end
    if p[i]>0
      j=1
      cump = d[i]
      k=p[i]
      active[i]=false
      while active[k] > 0
        j+=1
        cump *= d[k]
        active[k] = false
        k=p[k]
      end
      append!(eigv,[cump^(1.0/j)*exp(2*im*π*m/j) for m=1:j])
    else
      push!(eigv,0.0 + 0.0im)
    end
  end
  return eigv
end

pdeig(m::Matrix) = pdeig(pdfact(m)...)

n = 4   # testing vector to matrix transformation of permutations
σ=randperm(n)
@assert mat2perm(permmatrix(σ))==σ

例如：

m = [ 0.0 1.0 0.0 ; 2.0 0.0 0.0 ; 0.0 0.0 0.0 ]
pdeig(m)

输出：

    3-element Array{Complex{Float64},1}:
 -1.41421+1.73191e-16im 
  1.41421-3.46382e-16im 
              0.0+0.0im

由于这些矩阵是可对角化的，特征值应该提供对角矩阵（只需在它们上使用diagm）。

这些矩阵非常结构化，适当的 Julia 处理将为这些矩阵定义一个类型，然后定义各种线性代数函数以在该类型上调度。

如果出现错误，只需添加评论，我会尝试修复它们（或者如果我碰巧看到一个不错的重构，那么我会编辑）。

顺便说一句，计算会引入小的数值误差，这些应该不是问题，可以通过适当的舍入来消除（所以不必害怕 -1.0 是 -1.0+1.234234e-16im）