我正在尝试编写一个程序，该程序构造矩阵并对其执行奇异值分解。我正在评估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])
但我总是会收到广播错误，但不确定如何解决。
首先，在MATLAB中您要分配给D(:,i)，而在python中您要分配给D[i]。后者等效D[i, ...]于您的情况D[i, :]。相反，您似乎需要D[:, i]。
其次，在MATLAB中使用线性索引进入2d数组（即a和b）将为您提供扁平化的视图。如果使用numpy进行操作，则会得到数组的切片，就像我在提到的那样D[i]。
您可以通过附加（或重塑）您的和数组来消除循环播放和获取所需的2d数组：.ravelab
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
这种工作方式是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)
这是我有时使用的技巧，但是如果您不熟悉该模式，则很难阅读。
旁注：最好使用示例数据，其中尺寸最终将具有不同的大小，以便您注意到是否最终要换位。
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句