将 numpy 一维数组用于线性代数的最佳方法
Best approach to use numpy 1d arrays for linear algebra
在使用 numpy 进行线性代数操作时，我遇到了很多“小”问题，因为 numpy 处理“向量”或一维数组的方式在我看来造成了不一致的行为。
我的问题是，我是否在使用 numpy 数组进行线性代数时犯了一些明显的错误，或者这是否就是它的工作原理并且没有其他明显的方法可以做到这一点？
例如，假设我想对两个向量执行单变量 OLS。
import numpy as np
from numpy import linalg as la
y = np.arange(10)
x = np.arange(10)
print(x.shape)
ols = la.inv(x.T@x)@(x.T@y)
LinAlgError: 0-dimensional array given. Array must be at least two-dimensional
所以一种解决方案是强制数组具有额外的维度：
import numpy as np
from numpy import linalg as la
y = np.arange(10).reshape(-1, 1)
x = np.arange(10).reshape(-1, 1)
ols = la.inv(x.T@x)@(x.T@y)
print(ols)
>>> [[1.]]
然后以为问题解决了！但不完全是。如果我在 X 轴上有更多的维度，计算 t 值就成了问题。
y = np.arange(10).reshape(-1, 1)
X = np.arange(20).reshape(10, 2)
b_hat = la.inv((X.T@X))@(X.T@y)
# Calculate standard errors
residual = y - X@b_hat
sigma_hat = residual.T@residual/(y.size - b_hat.size)
b_var = sigma_hat*la.inv(X.T@X)
b_std = np.sqrt(b_var.diagonal()) # The diagonal method returns 1d array.
# Calculate t-values
t_values = b_hat/b_std
print(t_values)
>>> [[2.47854011e+13 2.67712930e+13]
[1.40888694e+00 1.52177182e+00]]
这当然不是故意的。为什么会这样？这是因为 np.sqrt(b_var.diagonal()) returns b_std 的 (2,) 形状。因此，当我划分 b_hat/b_std numpy 检查它们是否具有相同的形状时，它们不是（b_hat 具有 (2, 1) 形状），并且 numpy 不会进行“真正的”划分，但是一些其他类型的划分。
这个问题的解决方案当然是再次使用.reshape(-1, 1)，但是我的计算会越来越复杂，所以总是检查一个向量是作为一维数组还是二维数组返回是很麻烦的，并且然后重塑它。如果我不小心将矩阵重塑为向量，这也容易出错。
那么，我是不是在使用 numpy 数组进行线性代数方面犯了一些明显的错误，或者如果这就是它的工作原理并且没有其他明显的方法可以做到这一点？
也许记住这些规则可能有助于在使用 numpy 时提供您正在寻找的警告（术语轴和维度在这里含义相同）：
- 在数学中，当我们写下一系列数字
[1 2 3 4] 时，我们选择与该符号相关联的语义会随着上下文的不同而有所不同。有时我们将其视为单轴数组（这是正确的语义），但有时我们将其视为“1 行，4 列”。您还能如何证明数学家声称列向量在转置时给出行向量，反之亦然？术语“转置”表示行和列的交换，这本身意味着是两个轴，而不仅仅是一个。在 numpy 的情况下，同一事物的语义将 一致且严格地 是“长度为 4 的单轴”而不是“长度为 1 的第一轴和第二轴长度 4".
- 在
numpy 中，与数学的情况一样，转置的想法只有在至少有两个轴的情况下才有意义。如上所述，在数学中，我们没有一致的符号来区分单轴数组和双轴数组，因此这条规则确实没有实际意义。在 numpy 中执行 arr.T 只需 returns arr，如果 arr 恰好是一个单轴数组。
- 在
numpy中，我们可以在我们选择的任何位置添加一个额外的单位长度轴。一种表示法是 arr.reshape(n1,n2,...1,...,nk)（即在现有的以逗号分隔的轴长中间插入 1）。另一种方法是使用索引符号 arr[:,:,...,None,...,:]（即，使用与轴一样多的逗号分隔冒号，并在其中插入 None）。代替None，也可以用np.newaxis，只是比较冗长
- 基于以上内容，如果
arr 具有单轴（例如, 形状 (3,))。 （如何为单轴数组定义矩阵乘法？）相反，表达式 return 是 arr 和 arr.T 以及 return 的元素的和积将它作为一个标量（甚至 return 它不是一个单元素数组）。
- 在
numpy 中，当在两个相同形状的数组之间使用算术和比较运算符时，将“按元素”应用（这意味着在属于两个数组的每对对应元素之间） .这将产生一个新数组，其形状与操作数数组的形状相同。
- 对于算术运算符和比较运算符，如果两个操作数是形状不同但可广播为共同形状的数组，同样，运算符在广播后按元素应用，结果将再次是一个数组广播生成的形状。
- 对于算术运算符和比较运算符，如果其中一个操作数是标量，则该标量将被视为形状为
(1,) 的数组，然后将应用先前的（广播）规则。
- 虽然上述第 5、6、7 点实际上增加了
numpy 的表达能力，但它们经常给新学习者带来惊喜。例如，1.0 / arr 其中 arr 是 [1 2 3 4] 将生成一个由值 [1.0/1 1.0/2 1.0/3 1.0/4] 组成的新数组。 （我认为这是你进行除法时遇到的惊喜之一）
- 如果
arr具有(3,4,1,5,2,1,1)的形状，那么arr.squeeze()将去掉单位长度轴，因此return形成一个形状数组(3,4,5,2)
- 当我们索引一个多维数组时，我们通常期望结果具有较低的维度（更少的轴，and/or 相同轴的长度更小）或与被索引的数组相同的维度。在
numpy中，索引表达式如arr[my_index_arr]可以产生比原始数组arr更复杂、维数更高的形状。同样，这是一个强大的表达功能，有时可以 surprise/confuse 新学习者。在numpy中，这叫做Advanced Indexing with Integer Arrays
为了强调上面的一点，当你的阵列只有一个轴（形状像 (L,)）时，请特别注意你的期望。
即使在面向矩阵的 MATLAB 中（最初只有二维矩阵），我发现跟踪维度占调试工作的 80%。在 numpy 中，密切关注维度没有捷径可走。不要假设。
在你的第一种情况下，数组是 1d:
In [305]: x
Out[305]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [306]: x.T # same, only one axis to 'switch'
Out[306]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
一维数组的矩阵乘积是 dot/inner 乘积，结果是标量
In [307]: x.T@x
Out[307]: 285
@/matmul 不喜欢使用标量。 np.dot 对它们没问题，但是 matmul 是为处理 'batches' 数组而创建的：
In [308]: (x@x)@(x@y)
Traceback (most recent call last):
File "<ipython-input-308-2d486b202d47>", line 1, in <module>
(x@x)@(x@y)
TypeError: unsupported operand type(s) for @: 'numpy.int64' and 'numpy.int64'
然后创建二维数组，(10,1) 形状
In [309]: y = y[:,None]
In [310]: y
Out[310]:
array([[0],
[1],
...
[9]])
In [311]: y.T # (1,10) shape
Out[311]: array([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
现在 matmul (1,10) 与 (10,1) 在 10 上求和生成 (1,1) 数组：
In [312]: y.T@y
Out[312]: array([[285]])
In [313]: _.shape
Out[313]: (1, 1)
两个 (1,1) 数组与 @ 一起工作以再次生成 (1,1)，并且 inv 可以。
随着X
In [314]: X.shape
Out[314]: (10, 2)
In [315]: X.T@X # (2,10) with (10,2) producing (2,2)
Out[315]:
array([[1140, 1230],
[1230, 1330]])
In [316]: X.T@y # (2,10) with (10,1) producing (2,1)
Out[316]:
array([[570],
[615]])
In [317]: (X.T@X)@(X.T@y) # (2,2) with (2,1) producing (2,1)
Out[317]:
array([[1406250],
[1519050]])
所有这些的关键是 (K,n) 和 (n,M) 产生 (K,M) 以及共享 n 上的乘积总和 - 第一个的最后一个与第二个到最后一个第二个参数。那个高中跨行下列做矩阵乘积的方法。
matmul 文档讨论将 1d 数组提升为 2d 以执行相同的操作 - 但结果中删除了额外的维度。
下一个@将 (10,2) 与 (2,1) 连接起来生成 (10,1)（对 2 求和）：
In [319]: X@Out[317]
Out[319]:
array([[ 1519050],
[ 7369650],
[13220250],
...
[54174450]])
可以用(10,1)相减y.
residual.T@residual 是 (1,10) 与 (10,1) producint (1,1) (magnitude).
sigma_hat*la.inv(X.T@X) 是 (1,1) 乘以 (2,2) 得到 (2,2)（如果 sigma_hat 是标量，那效果也一样好。inv(A) 在 (1,1) 上只是 1/A.
b_var.diagonal() throws away 的值为 b_var。由于它是标量时间 X.T@X，我们可以检查：
In [323]: X.T@X
Out[323]:
array([[1140, 1230],
[1230, 1330]])
In [324]: (X.T@X).diagonal()
Out[324]: array([1140, 1330])
这与在大小 10 维度上求和相同
In [325]: (X*X).sum(0)
Out[325]: array([1140, 1330])
In [326]: np.einsum('ij,ij->j',X,X)
Out[326]: array([1140, 1330])
它将二维视为 'batch':
In [328]: [X[:,i]@X[:,i] for i in range(2)]
Out[328]: [1140, 1330]
matmul 旨在与 'batches' 一起使用 - 对于 3d（和更高）数组：
In [329]: Xt = X.T
In [330]: Xt.shape
Out[330]: (2, 10)
In [331]: Xt[:,None,:]@Xt[:,:,None] # (2,1,10) with (2,10,1)=>(2,1,1)
Out[331]:
array([[[1140]],
[[1330]]])
In [332]: (Xt[:,None,:]@Xt[:,:,None]).squeeze()
Out[332]: array([1140, 1330])
这表明 X 从一开始就应该是 (2,10) 甚至 (2,1,10)。 （和 y (1,10)？）。
这个答案有点啰嗦，但希望详细信息能让您了解如何跟踪维度。它旨在补充其他答案的一般原则。
在使用 numpy 进行线性代数操作时，我遇到了很多“小”问题，因为 numpy 处理“向量”或一维数组的方式在我看来造成了不一致的行为。
我的问题是，我是否在使用 numpy 数组进行线性代数时犯了一些明显的错误，或者这是否就是它的工作原理并且没有其他明显的方法可以做到这一点？
例如，假设我想对两个向量执行单变量 OLS。
import numpy as np
from numpy import linalg as la
y = np.arange(10)
x = np.arange(10)
print(x.shape)
ols = la.inv(x.T@x)@(x.T@y)
LinAlgError: 0-dimensional array given. Array must be at least two-dimensional
所以一种解决方案是强制数组具有额外的维度：
import numpy as np
from numpy import linalg as la
y = np.arange(10).reshape(-1, 1)
x = np.arange(10).reshape(-1, 1)
ols = la.inv(x.T@x)@(x.T@y)
print(ols)
>>> [[1.]]
然后以为问题解决了！但不完全是。如果我在 X 轴上有更多的维度，计算 t 值就成了问题。
y = np.arange(10).reshape(-1, 1)
X = np.arange(20).reshape(10, 2)
b_hat = la.inv((X.T@X))@(X.T@y)
# Calculate standard errors
residual = y - X@b_hat
sigma_hat = residual.T@residual/(y.size - b_hat.size)
b_var = sigma_hat*la.inv(X.T@X)
b_std = np.sqrt(b_var.diagonal()) # The diagonal method returns 1d array.
# Calculate t-values
t_values = b_hat/b_std
print(t_values)
>>> [[2.47854011e+13 2.67712930e+13]
[1.40888694e+00 1.52177182e+00]]
这当然不是故意的。为什么会这样？这是因为 np.sqrt(b_var.diagonal()) returns b_std 的 (2,) 形状。因此，当我划分 b_hat/b_std numpy 检查它们是否具有相同的形状时，它们不是（b_hat 具有 (2, 1) 形状），并且 numpy 不会进行“真正的”划分，但是一些其他类型的划分。
这个问题的解决方案当然是再次使用.reshape(-1, 1)，但是我的计算会越来越复杂，所以总是检查一个向量是作为一维数组还是二维数组返回是很麻烦的，并且然后重塑它。如果我不小心将矩阵重塑为向量，这也容易出错。
那么，我是不是在使用 numpy 数组进行线性代数方面犯了一些明显的错误，或者如果这就是它的工作原理并且没有其他明显的方法可以做到这一点？
也许记住这些规则可能有助于在使用 numpy 时提供您正在寻找的警告（术语轴和维度在这里含义相同）：
- 在数学中，当我们写下一系列数字
[1 2 3 4]时，我们选择与该符号相关联的语义会随着上下文的不同而有所不同。有时我们将其视为单轴数组（这是正确的语义），但有时我们将其视为“1 行，4 列”。您还能如何证明数学家声称列向量在转置时给出行向量，反之亦然？术语“转置”表示行和列的交换，这本身意味着是两个轴，而不仅仅是一个。在numpy的情况下，同一事物的语义将 一致且严格地 是“长度为 4 的单轴”而不是“长度为 1 的第一轴和第二轴长度 4". - 在
numpy中，与数学的情况一样，转置的想法只有在至少有两个轴的情况下才有意义。如上所述，在数学中，我们没有一致的符号来区分单轴数组和双轴数组，因此这条规则确实没有实际意义。在numpy中执行arr.T只需 returnsarr，如果arr恰好是一个单轴数组。 - 在
numpy中，我们可以在我们选择的任何位置添加一个额外的单位长度轴。一种表示法是arr.reshape(n1,n2,...1,...,nk)（即在现有的以逗号分隔的轴长中间插入1）。另一种方法是使用索引符号arr[:,:,...,None,...,:]（即，使用与轴一样多的逗号分隔冒号，并在其中插入None）。代替None，也可以用np.newaxis，只是比较冗长 - 基于以上内容，如果
arr具有单轴（例如, 形状(3,))。 （如何为单轴数组定义矩阵乘法？）相反，表达式 return 是arr和arr.T以及 return 的元素的和积将它作为一个标量（甚至 return 它不是一个单元素数组）。 - 在
numpy中，当在两个相同形状的数组之间使用算术和比较运算符时，将“按元素”应用（这意味着在属于两个数组的每对对应元素之间） .这将产生一个新数组，其形状与操作数数组的形状相同。 - 对于算术运算符和比较运算符，如果两个操作数是形状不同但可广播为共同形状的数组，同样，运算符在广播后按元素应用，结果将再次是一个数组广播生成的形状。
- 对于算术运算符和比较运算符，如果其中一个操作数是标量，则该标量将被视为形状为
(1,)的数组，然后将应用先前的（广播）规则。 - 虽然上述第 5、6、7 点实际上增加了
numpy的表达能力，但它们经常给新学习者带来惊喜。例如，1.0 / arr其中arr是[1 2 3 4]将生成一个由值[1.0/1 1.0/2 1.0/3 1.0/4]组成的新数组。 （我认为这是你进行除法时遇到的惊喜之一） - 如果
arr具有(3,4,1,5,2,1,1)的形状，那么arr.squeeze()将去掉单位长度轴，因此return形成一个形状数组(3,4,5,2) - 当我们索引一个多维数组时，我们通常期望结果具有较低的维度（更少的轴，and/or 相同轴的长度更小）或与被索引的数组相同的维度。在
numpy中，索引表达式如arr[my_index_arr]可以产生比原始数组arr更复杂、维数更高的形状。同样，这是一个强大的表达功能，有时可以 surprise/confuse 新学习者。在numpy中，这叫做Advanced Indexing with Integer Arrays
为了强调上面的一点，当你的阵列只有一个轴（形状像 (L,)）时，请特别注意你的期望。
即使在面向矩阵的 MATLAB 中（最初只有二维矩阵），我发现跟踪维度占调试工作的 80%。在 numpy 中，密切关注维度没有捷径可走。不要假设。
在你的第一种情况下，数组是 1d:
In [305]: x
Out[305]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [306]: x.T # same, only one axis to 'switch'
Out[306]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
一维数组的矩阵乘积是 dot/inner 乘积，结果是标量
In [307]: x.T@x
Out[307]: 285
@/matmul 不喜欢使用标量。 np.dot 对它们没问题，但是 matmul 是为处理 'batches' 数组而创建的：
In [308]: (x@x)@(x@y)
Traceback (most recent call last):
File "<ipython-input-308-2d486b202d47>", line 1, in <module>
(x@x)@(x@y)
TypeError: unsupported operand type(s) for @: 'numpy.int64' and 'numpy.int64'
然后创建二维数组，(10,1) 形状
In [309]: y = y[:,None]
In [310]: y
Out[310]:
array([[0],
[1],
...
[9]])
In [311]: y.T # (1,10) shape
Out[311]: array([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
现在 matmul (1,10) 与 (10,1) 在 10 上求和生成 (1,1) 数组：
In [312]: y.T@y
Out[312]: array([[285]])
In [313]: _.shape
Out[313]: (1, 1)
两个 (1,1) 数组与 @ 一起工作以再次生成 (1,1)，并且 inv 可以。
随着X
In [314]: X.shape
Out[314]: (10, 2)
In [315]: X.T@X # (2,10) with (10,2) producing (2,2)
Out[315]:
array([[1140, 1230],
[1230, 1330]])
In [316]: X.T@y # (2,10) with (10,1) producing (2,1)
Out[316]:
array([[570],
[615]])
In [317]: (X.T@X)@(X.T@y) # (2,2) with (2,1) producing (2,1)
Out[317]:
array([[1406250],
[1519050]])
所有这些的关键是 (K,n) 和 (n,M) 产生 (K,M) 以及共享 n 上的乘积总和 - 第一个的最后一个与第二个到最后一个第二个参数。那个高中跨行下列做矩阵乘积的方法。
matmul 文档讨论将 1d 数组提升为 2d 以执行相同的操作 - 但结果中删除了额外的维度。
下一个@将 (10,2) 与 (2,1) 连接起来生成 (10,1)（对 2 求和）：
In [319]: X@Out[317]
Out[319]:
array([[ 1519050],
[ 7369650],
[13220250],
...
[54174450]])
可以用(10,1)相减y.
residual.T@residual 是 (1,10) 与 (10,1) producint (1,1) (magnitude).
sigma_hat*la.inv(X.T@X) 是 (1,1) 乘以 (2,2) 得到 (2,2)（如果 sigma_hat 是标量，那效果也一样好。inv(A) 在 (1,1) 上只是 1/A.
b_var.diagonal() throws away 的值为 b_var。由于它是标量时间 X.T@X，我们可以检查：
In [323]: X.T@X
Out[323]:
array([[1140, 1230],
[1230, 1330]])
In [324]: (X.T@X).diagonal()
Out[324]: array([1140, 1330])
这与在大小 10 维度上求和相同
In [325]: (X*X).sum(0)
Out[325]: array([1140, 1330])
In [326]: np.einsum('ij,ij->j',X,X)
Out[326]: array([1140, 1330])
它将二维视为 'batch':
In [328]: [X[:,i]@X[:,i] for i in range(2)]
Out[328]: [1140, 1330]
matmul 旨在与 'batches' 一起使用 - 对于 3d（和更高）数组：
In [329]: Xt = X.T
In [330]: Xt.shape
Out[330]: (2, 10)
In [331]: Xt[:,None,:]@Xt[:,:,None] # (2,1,10) with (2,10,1)=>(2,1,1)
Out[331]:
array([[[1140]],
[[1330]]])
In [332]: (Xt[:,None,:]@Xt[:,:,None]).squeeze()
Out[332]: array([1140, 1330])
这表明 X 从一开始就应该是 (2,10) 甚至 (2,1,10)。 （和 y (1,10)？）。
这个答案有点啰嗦，但希望详细信息能让您了解如何跟踪维度。它旨在补充其他答案的一般原则。