Python实现MATLAB矩阵运算触发TypeError列表索引错误求助

报错原因 #

• 你定义的sx、sy、F为Python原生列表，不支持[i,k]这种二维元组索引，该语法是MATLAB的索引习惯，不适用于Python列表
• 三个列表未预先分配空间，直接按索引赋值会触发越界错误
• Python中*对列表是重复操作，对numpy数组是逐元素乘法，你需要的矩阵乘法要使用@运算符或np.dot()方法
• 原MATLAB逻辑中sx、sy是每次循环生成的临时向量（分别为1×4行向量、4×1列向量），不需要全局存储所有循环的sx/sy值

修正后可运行代码（完全对齐MATLAB逻辑） #

import numpy as np                                                                              

B = np.array([[-6.08066634428988e-10, -8.61023850910464e-11, 5.48222828615260e-12, -9.49229025004441e-14], 
    [-3.38148313553674e-11, 6.47759097087283e-12, 1.14900158474371e-13, -5.70078947874486e-15], 
    [-2.55893304237669e-13, -1.40941560399352e-13, 5.76510238931847e-15, -5.52980385181738e-17], 
    [3.39795122177475e-15, 7.95704191204353e-16, -5.31260642039813e-17, 7.83532802015832e-19]])                                                                        

# 生成网格，和MATLAB逻辑完全一致
X, Y = np.meshgrid(np.arange(0, 3, 0.01*3), np.arange(0, 15, 0.01*15))                           
# 预先分配结果矩阵空间，和X维度相同
F = np.zeros_like(X)
poly_order = B.shape[0]

for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        # 生成1×4的sx行向量
        sx = np.array([X[i,j]**k for k in range(poly_order)]).reshape(1, -1)
        # 生成4×1的sy列向量
        sy = np.array([Y[i,j]**l for l in range(poly_order)]).reshape(-1, 1)
        # 矩阵乘法得到标量结果赋值
        F[i,j] = sx @ B @ sy

向量化优化版本（去除双重循环，运行效率更高） #

import numpy as np                                                                              

B = np.array([[-6.08066634428988e-10, -8.61023850910464e-11, 5.48222828615260e-12, -9.49229025004441e-14], 
    [-3.38148313553674e-11, 6.47759097087283e-12, 1.14900158474371e-13, -5.70078947874486e-15], 
    [-2.55893304237669e-13, -1.40941560399352e-13, 5.76510238931847e-15, -5.52980385181738e-17], 
    [3.39795122177475e-15, 7.95704191204353e-16, -5.31260642039813e-17, 7.83532802015832e-19]])                                                                        

X, Y = np.meshgrid(np.arange(0, 3, 0.01*3), np.arange(0, 15, 0.01*15))
poly_order = B.shape[0]

# 批量生成所有位置的多项式特征
X_pows = np.stack([X**k for k in range(poly_order)], axis=-1)
Y_pows = np.stack([Y**l for l in range(poly_order)], axis=-1)[..., np.newaxis]

# 批量矩阵计算得到最终结果
F = (X_pows[..., np.newaxis, :] @ B @ Y_pows).squeeze()

内容的提问来源于stack exchange，提问作者Jacobo22