编写函数来处理向量和矩阵答案

【问题标题】：Writing functions to work on vectors and matrices编写函数来处理向量和矩阵
【发布时间】：2012-01-22 08:35:05
【问题描述】：

作为更大函数的一部分，我正在编写一些代码来生成一个向量/矩阵（取决于输入），其中包含输入向量/矩阵“x”的每一列的平均值。这些值存储在与输入向量形状相同的向量/矩阵中。

我的初步解决方案是在一维和矩阵数组上工作的非常（！）混乱：

# 'x' is of type array and can be a vector or matrix.
import scipy as sp
shp = sp.shape(x)
x_mean = sp.array(sp.zeros(sp.shape(x)))

try: # if input is a matrix
    shp_range = range(shp[1])
    for d in shp_range:
        x_mean[:,d] = sp.mean(x[:,d])*sp.ones(sp.shape(z))
except IndexError: # error occurs if the input is a vector
    z = sp.zeros((shp[0],))
    x_mean = sp.mean(x)*sp.ones(sp.shape(z))

来自 MATLAB 背景，这就是它在 MATLAB 中的样子：

[R,C] = size(x);
for d = 1:C,
    xmean(:,d) = zeros(R,1) + mean(x(:,d));
end

这适用于向量和矩阵，没有错误。

我的问题是，如何让我的 python 代码在没有（丑陋的）try/except 块的情况下同时处理向量和矩阵格式的输入？

谢谢！

【问题讨论】：

可能的欺骗：stackoverflow.com/q/8287047/48837
公平地说，matlab 根本没有一维数组。 matlab中的“一维”向量是二维的。
如果您打算使用 2D/1D 输出来计算标准，那么您不需要 2D 输出，并且可以通过广播在任何情况下使用 1D 数组。

标签： python matlab numpy scipy

【解决方案1】：

对于平均值计算本身，您不需要区分向量和矩阵 - 如果您使用 axis 参数 Numpy 将沿着向量（对于向量）或列（对于矩阵）执行计算。然后构建输出，你可以使用一个很好的老式列表理解，虽然它对于巨大的矩阵可能有点慢：

import numpy as np
m = np.mean(x,axis=0) # For vector x, calculate the mean. For matrix x, calculate the means of the columns
x_mean = np.array([m for k in x]) # replace elements for vectors or rows for matrices

使用列表推导式创建输出很慢，因为它必须分配内存两次——一次用于列表，一次用于数组。使用np.repeat 或np.tile 会更快，但对于向量输入来说会很有趣——输出将是一个嵌套矩阵，每行有一个 1 长的向量。如果速度比优雅更重要，您可以在以下情况下替换最后一行：

if len(x.shape) == 1:
    x_mean = m*np.ones(len(x))
else:
    x_mean = np.tile(m, (x.shape[1],1))

顺便说一句，您的 Matlab 代码对于行向量和列向量的行为不同（尝试使用 x 和 x' 运行它）。

【讨论】：

好点。我匆匆忙忙，应该花更多时间来解释事情。你的回答很好地解释了它。（实际上是正确的，不像我的:)）
@JoeKington - 谢谢！我根据您对速度的原始答案添加了一个替代方案。
对 MATLAB 代码中的差异感到抱歉 - 我只是想说明列向量的意义（我不希望我的程序中有任何行向量）。不过感谢您的澄清！

【解决方案2】：

首先关于 numpy 广播的快速说明。当我从 matlab 切换到 python 时，广播让我有点困惑，但是一旦我花时间去理解它，我意识到它是多么有用。要了解有关广播的更多信息，请查看http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html，

因为在 numpy 中广播 (m,) 数组（您称之为向量）本质上等同于 (1, m) 数组或 (1, 1, m) 数组等。似乎您希望 (m,) 数组的行为类似于 (m, 1) 数组。我相信有时会发生这种情况，尤其是在 linalg 模块中，但如果你要这样做，你应该知道你正在打破 numpy 约定。

有这个警告有代码：

import scipy as sp
def my_mean(x):
    if x.ndim == 1:
        x = x[:, sp.newaxis]
    m = sp.empty(x.shape)
    m[:] = x.mean(0)
    return sp.squeeze(m)

还有一个例子：

In [6]: x = sp.arange(30).reshape(5,6)

In [7]: x
Out[7]: 
array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11],
       [12, 13, 14, 15, 16, 17],
       [18, 19, 20, 21, 22, 23],
       [24, 25, 26, 27, 28, 29]])

In [8]: my_mean(x)
Out[8]: 
array([[ 12.,  13.,  14.,  15.,  16.,  17.],
       [ 12.,  13.,  14.,  15.,  16.,  17.],
       [ 12.,  13.,  14.,  15.,  16.,  17.],
       [ 12.,  13.,  14.,  15.,  16.,  17.],
       [ 12.,  13.,  14.,  15.,  16.,  17.]])

In [9]: my_mean(x[0])
Out[9]: array([ 2.5,  2.5,  2.5,  2.5,  2.5,  2.5])

这比使用tile快，时间如下：

In [1]: import scipy as sp

In [2]: x = sp.arange(30).reshape(5,6)

In [3]: m = x.mean(0)

In [5]: timeit m_2d = sp.empty(x.shape); m_2d[:] = m
100000 loops, best of 3: 2.58 us per loop

In [6]: timeit m_2d = sp.tile(m, (len(x), 1))
100000 loops, best of 3: 13.3 us per loop

【讨论】：