R语言是否有更快的方法生成符合指定规则的矩阵？

需求说明 #

生成n阶矩阵，规则如下（基于R的1索引规则）：

• 第i行的前min(i+1, n)个元素值为1/(i+1)
• 其余元素值为0
• 支持n最大为10000的大尺寸场景

最优实现方案 #

极简向量化实现（推荐，性能最高） #

利用R内置向量化操作，全程无循环，n=10000时运行耗时<1s：

n <- 10000 # 可替换为任意需要的尺寸
# 先构造每行统一取值的矩阵，按列填充自动让整行值相同
m <- matrix(1/(1:n + 1), nrow = n, ncol = n)
# 仅将列号 >= 行号+2的位置设为0，保留需要的非零区域
m[upper.tri(m, diag = -1)] <- 0

结果验证（以n=4为例） #

n <- 4
m <- matrix(1/(1:n + 1), nrow = n, ncol = n)
m[upper.tri(m, diag = -1)] <- 0
MASS::fractions(m)

输出结果完全符合要求：

[,1] [,2] [,3] [,4]
[1,] 1/2  1/2  0    0   
[2,] 1/3  1/3  1/3  0   
[3,] 1/4  1/4  1/4  1/4 
[4,] 1/5  1/5  1/5  1/5 

大尺寸内存优化方案（n>10000时推荐） #

如果需要进一步降低内存占用，可采用稀疏矩阵存储，内存占用仅为稠密矩阵的50%左右：

library(Matrix)
n <- 10000
# 计算每行非零元素个数
row_len <- pmin(1:n + 1, n)
# 构造非零元素的行、列索引与对应值
rows <- rep(1:n, times = row_len)
cols <- unlist(lapply(row_len, seq_len))
vals <- 1/(rows + 1)
# 生成稀疏矩阵
m <- sparseMatrix(i = rows, j = cols, x = vals, dims = c(n, n))

原代码问题说明 #

你之前使用upper.tri(m, diag=FALSE)会将所有列号>行号的位置设为0，丢失了列号=行号+1的次对角线元素。将diag参数设置为-1后，upper.tri仅匹配列号>=行号+2的位置，即可保留需要的非零区域。

内容的提问来源于stack exchange，提问作者Homer Jay Simpson