嘿，这个问题问到点子上了！循环实现虽然直观易懂，但在处理大规模数据时，Python层面的循环会拖慢效率。NumPy的核心优势就是向量化运算——用底层C实现的操作替代Python循环，速度能提升好几个量级。这里给你两种最快的实现方式：

方法一：用np.einsum（灵活且高效） #

np.einsum通过爱因斯坦求和约定，能非常灵活地处理多维数组的运算，完美适配这个场景：

import numpy as np
diagonals = np.array([[1,2,3],[4,5,6],[7,8,9]])
# 一步生成三维对角矩阵数组
diag_matrices = np.einsum('...i,ij->...ij', diagonals, np.eye(diagonals.shape[1]))
print(diag_matrices)

简单解释下：...i匹配diagonals的任意前置维度+最后一维，ij是单位矩阵的维度，通过求和约定，每个对角线行都会和单位矩阵对应相乘，自动扩展成三维数组，完全没有循环开销。

方法二：维度扩展+广播乘法（更直观） #

如果你对einsum的语法不太熟悉，用维度扩展结合广播的方式也能达到同样的效果：

import numpy as np
diagonals = np.array([[1,2,3],[4,5,6],[7,8,9]])
n = diagonals.shape[1]
# 把diagonals从(3,3)扩展为(3,3,1)，再和(3,3)的单位矩阵广播相乘
diag_matrices = np.expand_dims(diagonals, axis=-1) * np.eye(n)
print(diag_matrices)

原理很简单：维度扩展后，每个对角线元素会沿着新增的维度“复制”，然后和单位矩阵相乘时，只有对角位置的元素会保留（单位矩阵的对角是1，其他是0），刚好生成我们需要的对角矩阵。

为什么这两种方法更好？ #

当你的diagonals规模变大（比如有1000行，每行1000个元素），这两种向量化方法的速度会比Python循环快几十甚至上百倍——因为NumPy把运算放到了底层C代码中执行，避开了Python循环的解释开销。而且代码更简洁，可读性也更强。

你可以用np.array_equal验证一下，这两种方法生成的结果和你原来的循环实现完全一致：

# 循环实现的结果
diag_matrices_loop = np.zeros([3,3,3])
for i in range(3):
    diag_matrices_loop[i] = np.diag(diagonals[i,:])

print(np.array_equal(diag_matrices, diag_matrices_loop))  # 输出True

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