时间戳

首页 > 解决方案 > 如何区分 Python 或 Matlab 是否错误/故障？

python - 如何区分 Python 或 Matlab 是否错误/故障？

问题描述

我正在尝试使用 SVD 和特征分解进行一些使用动态模式分解的数据分析。我遇到了一个从 Matlab 和 Python 获得不同结果的简单问题。我很困惑，不知道为什么 Python 给了我错误的结果/矩阵值，但一切看起来（我认为是）都是正确的。

所以这次我没有使用真实数据并查看大型数据集，而是生成了数据。我将尝试查看特征分解后的特征值图。我还对数据使用了延迟嵌入，因为我将使用只有 (2x100) 的数据向量，因此我将执行一种 Hankel 矩阵来丰富具有 10 个延迟的数据。

clear all; close all; clc;
data = linspace(1,100);
data2 = linspace(2,101);
data = [data;data2];
numDelays = 10;
relTol= 10^-6;

%% Create first and second snap shot matrices for DMD. Any columns with missing
% data are not used.
disp('Constructing Data Matricies:')
X = zeros((numDelays+1)*size(data,1),size(data,2)-(numDelays+1));
Y = zeros(size(X));

for i = 1:numDelays+1
   X(1 + (i-1)*size(data,1):i*size(data,1),:) = ...
       data(:,(i):size(data,2)-(numDelays+1) + (i-1));
   Y(1 + (i-1)*size(data,1):i*size(data,1),:) = ...
       data(:,(i+1):size(data,2)-(numDelays+1) + (i));
end
[U,S,V] = svd(X);
r = find(diag(S)>S(1,1)*relTol,1,'last');
disp(['DMD subspace dimension:',num2str(r)])
U = U(:,1:r);
S = S(1:r,1:r);
V = V(:,1:r);
Atil = (U'*Y)*V*(S^-1);
[what,lambda] = eig(Atil);
Phi = (Y*V)*(S^-1)*what;

Keigs = diag(lambda);
tt = linspace(0,2*pi,101);
figure;
plot(real(Keigs),imag(Keigs),'ro')
hold on
plot(cos(tt),sin(tt),'--')

import scipy.io as sc
import math as m
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sys
from numpy import dot, multiply, diag, power, pi, exp, sin, cos, cosh, tanh, real, imag
from scipy.linalg import expm, sinm, cosm, fractional_matrix_power, svd, eig, inv


def dmd(X, Y, relTol):
    U2,Sig2,Vh2 = svd(X, False) # SVD of input matrix
    S = np.zeros((Sig2.shape[0], Sig2.shape[0]))  # Create S matrix with zeros based on Diag of S
    np.fill_diagonal(S, Sig2)  # Fill diagonal of S matrix with the nonzero values
    r = np.count_nonzero(np.diag(S) > S[0,0] * relTol) # rank truncation
    U = U2[:,:r]
    Sig = diag(Sig2)[:r,:r]     #GOOD =)
    V = Vh2.conj().T[:,:r]
    Atil = dot(dot(dot(U.conj().T, Y), V), inv(Sig)) # build A tilde
    print(Atil)
    mu,W = eig(Atil)
    Phi = dot(dot(dot(Y, V), inv(Sig)), W) # build DMD modes
    return mu, Phi

data = np.array([(np.linspace(1,100,100)),(np.linspace(2,101,100))])
Data = np.array(data)
#########           Choose number of Delays         ###########
# observable (coordinates of feature points). Setting to zero means only
# experimental observables will be used.
numDelays = 10
relTol = 10**-6
##########      Create Data Matrices for DMD        ###############
# Create first and second snap shot matrices for DMD. Any columns with missing
# data are not used.
X = np.zeros(((numDelays + 1) * data.shape[0], data.shape[1] - (numDelays + 1)))
Y = np.zeros(X.shape)

for i in range(1, numDelays + 2):
    X[0 + (i - 1) * Data.shape[0]:i * Data.shape[0], :] = Data[:, (i):Data.shape[1] - (numDelays + 1) + (i - 0)]

    Y[0 + (i - 1) * Data.shape[0]:i * Data.shape[0], :] = Data[:, (i + 0):Data.shape[1] - (numDelays + 1) + (i)]
Keigs, Phi = dmd(X, Y, relTol)

tt = np.linspace(0,2*np.pi,101)
plt.figure()
plt.plot(np.cos(tt),np.sin(tt),'--')
plt.plot(Keigs.real,Keigs.imag,'ro')
plt.title('DMD Eigenvalues')
plt.xlabel(r'Real $\ lambda$')
plt.ylabel(r'Imaginary $\ lambda$')
# plt.axes().set_aspect('equal')
plt.show()

所以在matlab和python中，我得到我的特征值都坐在单位圆上（正如预期的那样），我得到了一个，坐在1。

因此，当我查看 SVD 中的矩阵时，问题就来了，它们似乎具有不同的值。唯一相同的矩阵是“S 或 Sig”矩阵。其余的将不同的数字或 +/- 符号。最让我感兴趣的是 Atil 矩阵。在 matlab 中，它看起来像，[1.0157, -0.3116; 7.91229e-4, 0.9843] 和 python 它看起来像，[1.0, -4.508e-15; -4.439e-18, 1.0]

现在这可能由于数值错误而看起来略有偏差，但是当我查看真实数据并且这些不同时，它会打乱我的分析。

标签： pythonmatlabsvdeigenvalueeigenvector

解决方案

非方阵的 SVD 在 U 和 V 中不是唯一的。即使您有一个具有非零、非退化奇异值的方阵，U 和 V 中的奇异向量也仅在符号因子内是唯一的。 https://math.stackexchange.com/questions/644327/how-unique-on-non-unique-are-u-and-v-in-singular-value-decomposition-svd

此外，Matlab (LAPACK + BLAS) 和 scipy.linalg.svd 可能对 SVD 使用不同的算法。这可能会导致您所经历的差异。

推荐阅读