如何区分是Python还是Matlab是wrong/faulty?
How to distinguish if Python or Matlab is wrong/faulty?
我正在尝试使用 SVD 和特征分解来使用动态模式分解进行一些数据分析。我 运行 遇到了一个从 Matlab 和 Python 得到不同结果的简单问题。我很困惑,不知道为什么 Python 给我错误的 results/matrix 值,但一切看起来(我认为是)正确的。
因此,这次我没有使用真实数据并查看大型数据集,而是生成了数据。我将尝试查看特征分解后的特征值图。我还对数据使用了延迟嵌入,因为我将使用仅为 (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 or Sig' 矩阵。其余的将与数字或 +/- 符号不同。最让我感兴趣的是 Atil 矩阵。
在 matlab 中,它看起来像,
[1.0157, -0.3116; 7.91229e-4, 0.9843]
python 看起来,
[1.0, -4.508e-15; -4.439e-18, 1.0]
现在这可能由于数值错误而看起来略有偏差,但是当我查看实际数据时发现这些数据不同,这会扰乱我的分析。
non-square 矩阵的 SVD 在 U 和 V 中不是唯一的。即使您的方阵具有 non-zero、non-degenerate 奇异值,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 使用不同的算法。
这可能会导致您所经历的差异。
我正在尝试使用 SVD 和特征分解来使用动态模式分解进行一些数据分析。我 运行 遇到了一个从 Matlab 和 Python 得到不同结果的简单问题。我很困惑,不知道为什么 Python 给我错误的 results/matrix 值,但一切看起来(我认为是)正确的。
因此,这次我没有使用真实数据并查看大型数据集,而是生成了数据。我将尝试查看特征分解后的特征值图。我还对数据使用了延迟嵌入,因为我将使用仅为 (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 or Sig' 矩阵。其余的将与数字或 +/- 符号不同。最让我感兴趣的是 Atil 矩阵。 在 matlab 中,它看起来像, [1.0157, -0.3116; 7.91229e-4, 0.9843] python 看起来, [1.0, -4.508e-15; -4.439e-18, 1.0]
现在这可能由于数值错误而看起来略有偏差,但是当我查看实际数据时发现这些数据不同,这会扰乱我的分析。
non-square 矩阵的 SVD 在 U 和 V 中不是唯一的。即使您的方阵具有 non-zero、non-degenerate 奇异值,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 使用不同的算法。 这可能会导致您所经历的差异。