python - 如何将带有网格网格和数组的 Matlab 代码转换为 Python 代码？

Question

我正在尝试编写一个程序来构造一个矩阵并对其执行奇异值分解。我正在网格上评估函数 ax^2 +bx + 1。然后，我制作一个由 a 和 b 组成的统一网格。矩阵的行对应于不同的二次系数，而每列对应于计算函数的网格点。

matlab代码在这里:

% Collect data

x = linspace(-1,1,100);

[a,b] = meshgrid(0:0.1:1,0:0.1:1);

D=zeros(numel(x),numel(a));
sz = size(D)
% Build “Dose” matrix

for i=1:numel(a)

D(:,i) = a(i)*x.^2+b(i)*x+1;

end

% Do the SVD:

[U,S,V]=svd(D,'econ');

D_reconstructed = U*S*V';

plot(diag(S))
scatter3(a(:),b(:),V(:,1))

这是我尝试的解决方案:

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-1, 1, 100)

def f(x, a, b):
    return a*x*x + b*x + 1

a, b = np.mgrid[0:1:0.1,0:1:0.1]
#a = b = np.arange(0,1,0.01)

D = np.zeros((x.size, a.size))

for i in range(a.size):
    D[i] = a[i]*x*x +b[i]*x +1

U, S, V = np.linalg.svd(D)

plt.plot(np.diag(S))

fig = plt.figure()
ax = plt.axes(projection="3d")


ax.scatter(a, b, V[0])

但我总是遇到广播错误，我不知道如何修复。

Answer 1

首先，在 MATLAB 中，您要分配给 D(:,i)，但在 Python 中，您要分配给 D[i]。后者相当于 D[i, ...] ，在您的情况下是 D[i, :] 。相反，您似乎需要 D[:, i]。

其次，在 MATLAB 中使用二维数组的线性索引(即 a 和 b)将为您提供扁平 View 。如果您使用 numpy 执行此操作，您将得到数组的切片，正如我在 D[i] 中提到的那样。

您可以使用 broadcasting 消除循环并通过 .ravel 修改(或 reshape )您的 a 和 b 数组来获取所需的二维数组:

x = np.linspace(-1, 1, 100)[:, None]  # inject trailing singleton for broadcasting
a, b = np.mgrid[0:1:0.1, 0:1:0.1]
D = a.ravel() * x**2 + b.ravel() * x + 1

其工作方式是，在我们注入(inject)尾随单例之后，x 的形状为 (100, 1)(在 MATLAB 中，尾随单例是隐含的，在 numpy 中是前导单例) ，并且 a.ravel() 和 b.ravel() 都具有与 兼容的形状 (10*10,) (1, 10*10)，使广播成为可能 (100, 10*10) 形状。您还可以将对 ravel 的调用替换为

a, b = np.mgrid[...].reshape(2, -1)

这是我有时会使用的一个技巧，但如果您不熟悉该模式，则很难阅读。

旁注:最好使用尺寸最终大小不同的示例数据，以便您注意到某些内容是否最终被转置。