python - 从一个矩阵定义一个对角矩阵,其中包含子对角线的乘积
问题描述
我如何从矩阵 A 构造对角矩阵 D 其中:
- 对角矩阵 D 的第一个元素必须是矩阵 A 的下对角线中所有元素的乘积
- 对角矩阵 D 中的第二个元素必须是矩阵 A 的下对角线中除第一个元素之外的所有元素的乘积
- 对角矩阵 D 中的第三个元素必须是矩阵 A 的下对角线中除第一个和第二个元素之外的所有元素的乘积
- ...并且对角矩阵 D 的最后一个元素必须是 1
解决方案
这是一种方法:
def subdiag_prod_(a):
sub_diag = np.diagonal(a, offset=-1)
mask = np.triu(np.ones((sub_diag.shape*2))).astype('bool')
m = sub_diag[:,None].T * mask
ma = np.ma.array(sub_diag[:,None].T * mask, mask=~mask)
diag = np.prod(ma, axis=1).data
out = np.diag(diag)
last_row = np.zeros([out.shape[0]+1]*2)
last_row[:out.shape[0], :out.shape[1]] += out
return last_row
a = np.random.randint(1,5,(10,10))
array([[2, 2, 1, 4, 3, 1, 3, 1, 4, 4],
[2, 2, 2, 1, 1, 2, 3, 2, 2, 2],
[4, 2, 2, 4, 2, 1, 3, 3, 3, 4],
[2, 3, 3, 1, 1, 1, 4, 2, 3, 4],
[3, 3, 1, 1, 2, 1, 3, 4, 4, 3],
[1, 4, 4, 1, 1, 4, 1, 1, 1, 4],
[4, 2, 3, 2, 1, 4, 4, 1, 3, 2],
[4, 2, 2, 4, 4, 4, 1, 4, 3, 1],
[1, 3, 1, 1, 2, 2, 2, 3, 4, 1],
[1, 3, 2, 2, 3, 4, 1, 3, 2, 1]])
subdiag_prod(a)
array([[288., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 144., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 72., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 24., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 24., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 24., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 6., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 6., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 2., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]])
细节
第一步是采用 ndarray 的对角线np.diagonal
:
sub_diag = np.diagonal(a, offset=-1)
# array([2, 2, 3, 1, 1, 4, 1, 3, 2])
mask
我们可以通过 using创建一个np.tril
,然后我们可以使用它以指定的方式获取子对角元素的乘积:
mask = np.triu(np.ones((sub_diag.shape*2))).astype('bool')
现在我们可以ndarray
通过将掩码和次对角线相乘来创建一个掩码数组,使用上面的掩码:
mask = np.ma.array(sub_diag[:,None].T * mask, mask=~mask)
现在我们可以获取掩码数组的逐行乘积:
d = np.prod(ma, axis=1).data
# array([288, 144, 72, 24, 24, 24, 6, 6, 2])
然后简单地从中构建一个对角矩阵:
out = np.diag(d)
last_row = np.zeros([out.shape[0]+1]*2)
last_row[:out.shape[0], :out.shape[1]] += out
推荐阅读
- scala - 多项目路由在 playframework 上不起作用
- testing - 测试修改请求的golang中间件
- r - bwplot 中线而不是点
- reinforcement-learning - 强化学习,如何调整推荐系统
- java - 为命令行界面组织命令的好方法是什么?
- php - Doctrine querybuilder - 在原则中选择多对多连接的查询
- oozie - Oozie fork 多次调用相同的操作
- sql - 在 SQL Server 中将 RTF 转换为 varchar
- r - 将 2 个列表中的元素组合到 R 中的同一子集中
- python - 使用数值比较快速交错 numpy 数组