批量线性运算#
SciPy 的几乎所有线性代数函数现在都支持 N 维数组输入。这些操作尚未从数学上推广到高阶张量；相反，所指示的操作是针对输入标量、向量和/或矩阵的批次（或“堆栈”）执行的。
考虑 linalg.det 函数，它将矩阵映射到标量。
import numpy as np
from scipy import linalg
A = np.eye(3)
linalg.det(A)
np.float64(1.0)
有时我们需要相同维度的矩阵批次的行列式。
batch = [i*np.eye(3) for i in range(1, 4)]
batch
[array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]]),
array([[2., 0., 0.],
[0., 2., 0.],
[0., 0., 2.]]),
array([[3., 0., 0.],
[0., 3., 0.],
[0., 0., 3.]])]
我们可以通过循环或列表推导式对批次中的每个元素执行操作
[linalg.det(A) for A in batch]
[np.float64(1.0), np.float64(8.0), np.float64(27.0)]
然而，就像我们可能首先使用 NumPy 的广播和向量化规则来创建矩阵批次一样
i = np.arange(1, 4).reshape(-1, 1, 1)
batch = i * np.eye(3)
batch
array([[[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]],
[[2., 0., 0.],
[0., 2., 0.],
[0., 0., 2.]],
[[3., 0., 0.],
[0., 3., 0.],
[0., 0., 3.]]])
我们也可能希望在一次函数调用中对所有矩阵执行行列式运算。
linalg.det(batch)
array([ 1., 8., 27.])
在 SciPy 中，我们更喜欢使用“批次”（batch）而不是“堆栈”（stack），因为这个概念已推广到 N 维批次。假设输入是一个 2 x 4 的 3 x 3 矩阵批次。
batch_shape = (2, 4)
i = np.arange(np.prod(batch_shape)).reshape(*batch_shape, 1, 1)
input = i * np.eye(3)
在这种情况下，我们称批次形状为 (2, 4)，输入的核心形状为 (3, 3)。输入的总形状是批次形状和核心形状的总和（连接）。
input.shape
(2, 4, 3, 3)
由于每个 3 x 3 矩阵都被转换为零维标量，我们称输出的核心形状为 ()。输出的形状是批次形状和核心形状的总和，因此结果是一个 2 x 4 数组。
output = linalg.det(input)
output
array([[ 0., 1., 8., 27.],
[ 64., 125., 216., 343.]])
output.shape
(2, 4)
并非所有线性代数函数都映射到标量。例如，scipy.linalg.expm 函数将矩阵映射到具有相同形状的矩阵。
A = np.eye(3)
linalg.expm(A)
array([[2.71828183, 0. , 0. ],
[0. , 2.71828183, 0. ],
[0. , 0. , 2.71828183]])
在这种情况下，输出的核心形状是 (3, 3)，因此当批次形状为 (2, 4) 时，我们期望的输出形状是 (2, 4, 3, 3)。
output = linalg.expm(input)
output.shape
(2, 4, 3, 3)
将这些规则推广到具有多个输入和输出的函数也很简单。例如，scipy.linalg.eig 函数默认产生两个输出：一个向量和一个矩阵。
evals, evecs = linalg.eig(A)
evals.shape, evecs.shape
((3,), (3, 3))
在这种情况下，输出向量的核心形状是 (3,)，输出矩阵的核心形状是 (3, 3)。每个输出的形状仍然是批次形状加上核心形状。
evals, evecs = linalg.eig(input)
evals.shape, evecs.shape
((2, 4, 3), (2, 4, 3, 3))
当存在多个输入时，如果输入形状相同，则不会出现复杂情况。
evals, evecs = linalg.eig(input, b=input)
evals.shape, evecs.shape
((2, 4, 3), (2, 4, 3, 3))
当形状不相同时，规则也遵循逻辑。只要形状根据 NumPy 的广播规则 可广播，每个输入都可以有自己的批次形状。净批次形状是个别批次形状的广播结果，每个输出的形状是净批次形状加上其核心形状。
rng = np.random.default_rng(2859239482)
# Define input core shapes
m = 3
core_shape_a = (m, m)
core_shape_b = (m, m)
# Define broadcastable batch shapes
batch_shape_a = (2, 4)
batch_shape_b = (5, 1, 4)
# Define output core shapes
core_shape_evals = (m,)
core_shape_evecs = (m, m)
# Predict shapes of outputs: broadcast batch shapes,
# and append output core shapes
net_batch_shape = np.broadcast_shapes(batch_shape_a, batch_shape_b)
output_shape_evals = net_batch_shape + core_shape_evals
output_shape_evecs = net_batch_shape + core_shape_evecs
output_shape_evals, output_shape_evecs
((5, 2, 4, 3), (5, 2, 4, 3, 3))
# Check predictions
input_a = rng.random(batch_shape_a + core_shape_a)
input_b = rng.random(batch_shape_b + core_shape_b)
evals, evecs = linalg.eig(input_a, b=input_b)
evals.shape, evecs.shape
((5, 2, 4, 3), (5, 2, 4, 3, 3))
有一些函数，其参数或输出的核心维度（即核心形状的长度）可以是 1 或 2。在这些情况下，如果数组只有一个维度，则核心维度被视为 1；如果数组有两个或更多维度，则被视为 2。例如，考虑对 scipy.linalg.solve 的以下调用。最简单的情况是单个方阵 A 和单个向量 b
A = np.eye(5)
b = np.arange(5)
linalg.solve(A, b)
array([0., 1., 2., 3., 4.])
在这种情况下，A 的核心维度是 2（形状 (5, 5)），b 的核心维度是 1（形状 (5,)），输出的核心维度是 1（形状 (5,)）。
然而，b 也可以是一个二维数组，其中列被视为一维向量。
b = np.empty((5, 2))
b[:, 0] = np.arange(5)
b[:, 1] = np.arange(5, 10)
linalg.solve(A, b)
array([[0., 5.],
[1., 6.],
[2., 7.],
[3., 8.],
[4., 9.]])
b.shape
(5, 2)
乍一看，b 的核心形状似乎仍然是 (5,)，而我们只是以 (2,) 的批次形状执行了操作。然而，如果是这种情况，b 的批次形状将*前置*到核心形状，导致 b 和输出的形状为 (2, 5)。更仔细地思考，将两个输入和输出的核心维度都视为 2 是正确的；此时批次形状是 ()。
同样，每当 b 具有超过两个维度时，b 和输出的核心维度都被视为 2。例如，要解决批处理的三个完全独立的线性系统，每个系统只有一个右侧，b 必须作为三维数组提供：一个维度用于批次形状（(3,)），两个维度用于核心形状（(5, 1)）。
A = rng.random((3, 5, 5))
b = rng.random((3, 5, 1)) # batch shape (3,), core shape (5, 1)
linalg.solve(A, b).shape
(3, 5, 1)