julia - 对角化稀疏酉矩阵

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)