python - 跟踪1参数矩阵族的特征向量

Question

我的问题是这样的：我正在尝试通过（截断的）Karhunen-Loeve变换对随机过程进行频谱分解，但是我的协方差矩阵实际上是一个1参数矩阵族，我需要一种方法来估算/可视化我的随机过程取决于此参数。为此，我需要一种跟踪numpy.linalg.eigh（）产生的特征向量的方法。

为了让您了解我的问题，这里有一个示例玩具问题：假设我有一组点{xs}，并且随机过程R的协方差C（x，y）= 1 /（1 + a *（xy） ^ 2）取决于参数a。对于[0,1]范围内的网格点的随机样本和给定的a（例如a = 1）选择，我可以填充协方差矩阵并使用以下方法实现Karhunen-Loeve变换：

num_x = 1000
xs = numpy.array([numpy.random.uniform() for i in range(num_x)])
z=numpy.random.standard_normal(num_x)
a0=1
def cov(x,y,a=a0): return 1/(1+a*(x-y)^2)
cov_m = numpy.array(map(lambda y: map(lambda x: cov(x,y),xs),xs))
w,v=numpy.linalg.eigh(cov_m)
R=numpy.dot(v,z*w^0.5)

这将使我认识到R的值是在每个随机网格点xs处定义的。但是，我需要做的是-对于特定的实现（这意味着我的网格xs和我的随机系数z的特定选择）-跟踪R如何相对于我的协方差函数中的参数a变化。

如果我可以象征性地计算协方差矩阵并在事后指定a，这将很容易做到。但是，对于大型矩阵，这不是一个可行的选择。另一种方法是找到一种方法来跟踪numpy.linalg.eigh（）返回的每个特征向量。不幸的是，numpy似乎正在对它们进行重新排序，以便始终总是首先列出最小的特征值。这意味着，当我改变a时，特征向量将无法预测地重新排序，并且点积numpy.dot（v，z * w ^ 0.5）不再将相同的系数分配给相同的特征向量。

有没有解决的办法？

（这是来自ASKSAGE的交叉发布。我的实现使用的是Sage，但是由于问题不是Sage特定的，而且这里似乎有更多活动，所以我认为我会重新发布。如果不能接受交叉发布，我深表歉意；如果是这样，请删除它。）

编辑：基于下面的对话，我看到我需要添加有关此问题的性质的更多详细信息。

Karhunen-Loeve变换的思想是从光谱上分解随机过程R（x），如下所示：

R（x）= \ sum_ {i} Z_i \ lambda_i ^ {1/2} \ phi ^ {（i）}（x），

其中每个Z_i是一个i.d.具有标准正态分布的随机变量，每个\ phi ^ {{i}}（x）是x的非随机函数，由协方差矩阵的特征向量的解确定，每个\ lambda_i是与\ phi ^ {（i）}。

为了证明R对参数a的依赖性，我需要能够为每个\ phi ^ {（i）}唯一地分配一个系数Z_i。由于所有Z的分布都是相同的，所以该分配方式无关紧要，但是分配必须唯一，并且它不能依赖于\ lambda_i的对应值，因为这将取决于参数a。

解析解决方案很简单：只需为任意a计算特征向量和特征值，通过指定我的Z选择特定的R实现，然后观察R的实现如何依赖a。 （对于玩具模型，R通常会随着增加而迅速变化。）

这里的问题是如何以数字方式实现这一点。 Numpy显然在计算本征向量时会加扰，因此无法唯一地标记它们。我猜想这样做的唯一方法是深入研究底层代码，并专门实现某种类型的任意标记功能。

简而言之，我需要一种对numpy产生的特征向量进行排序的方法，该方法不依赖于相关特征值的大小。

有没有办法做到这一点？

更新：
我已经设法找到了这个问题的部分答案，但是我需要做更多的研究才能得出完整的答案。看来我将必须为此实现自己的算法。

对于我正在考虑的类型的矩阵，Lanczos算法（对于给定的初始向量选择）是确定性的，并且步骤不取决于我的参数选择。这给了我一个对称的三对角矩阵来求解特征值。

分而治之可能在这里起作用。似乎我可以实现它的一个版本，使我可以独立于相关的特征值来跟踪特征向量。至少可以以确定性，与参数无关的方式来实现“分隔”部分，但是我需要进一步了解该算法的“征服”部分才能确定。

Answer 1

经过一些研究，我设法为这个问题提出了两个部分答案。

首先是，对于没有零本征向量的实对称矩阵（可能还需要指定非简并矩阵），应该可行的是生成一种求解本征对问题的算法，该算法以固定顺序生成本征对，而与矩阵的选择。给定恒定的起始向量，Lanczos algorithm将以确定性方式为任意实数对称矩阵生成一个三对角矩阵。 divide-and-conquer algorithm的“ divide”部分具有类似的确定性，这意味着算法的唯一迭代次数取决于矩阵元素的值的部分是“ conquer”部分，从而求解了世俗方程：

1+ \ sum_j ^ m w_j ^ 2 /（d_j- \ lambda）= 0

因此，对于每个2x2块，问题归结为如何以一种实际上不依赖于原始矩阵值的方式对世俗方程的两个根进行排序。

第二部分解决方案更容易实现，但更容易失败。回想起来，这也是显而易见的。

同一矩阵的两个不同特征向量将始终彼此正交。因此，如果特征向量根据单个参数α平稳变化，则：

v_i（a）.v_j（a + da）= \ delta_ {ij} + O（da）

因此，随着参数α的变化，这给出了特征向量之间的自然映射。

这与David Zwicker和jorgeca提出的测量特征向量对之间的全局距离的想法类似，但是更容易实现。但是，在特征向量快速变化的区域或参数a的变化太大的区域中，此方法的实现将容易失败。

同样，在特征值相交处发生什么的问题也很有趣，因为在每个这样的相交处，系统都会退化。然而，在跨越简并性的允许特征向量集合中，将存在两个满足点积条件且可用作跨越简并性的基础，从而维持映射的两个向量。

当然，这是假设将特征向量视为参数空间上的光滑连续函数是正确的，但我不确定（正如jorgeca指出的那样）。