从 lower/upper 三角矩阵中给出一个元素
Giving an element from lower/upper triangular matrix
我有一个矩阵的上三角部分,主对角线存储为线性数组,如何从数组的线性索引中提取矩阵元素的 (i,j) 索引?
例如线性数组:[a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10]
是矩阵的存储
a0 a1 a2 a3
0 a4 a5 a6
0 0 a7 a8
0 0 0 a10
我找到了这个问题的解决方案,但没有主对角线:
index = (n*(n-1)/2) - (n-i)*((n-i)-1)/2 + j - i - 1
同样问题的解决方案,但对于具有对角线的下三角矩阵:
index = ((i + 1) * i / 2 + i).
此致,
我找到答案了!它是:
index = (n*(n+1)/2) - (n-i)*((n-i)+1)/2 + j - i
我的解法可能等同于,我没查过:
index = N * i - ((i - 1) * i) / 2 + (j - i)
这是一个完整的 Python 测试。我使用 Python 因为 Numpy 有 triu_indices
,它给出了上三角索引。
import numpy as np
def mksquare(N):
"""Make a square N by N matrix containing 0 .. N*N-1"""
return np.arange(N * N).reshape(N, N)
def mkinds(N):
"""Return all triu indexes for N by N matrix"""
return [(i,j) for i in range(N) for j in range(N) if i <= j]
def ij2linear(i, j, N):
"""Convert (i,j) 2D index to linear triu index for N by N array"""
return N * i - ((i - 1) * i) // 2 + (j - i)
def test(N):
"""Make sure my `mkinds` works for given N"""
arr = mksquare(N)
vec = arr[np.triu_indices(N)]
inds = mkinds(N)
expected = [arr[i, j] for (i, j) in inds]
actual = [vec[ij2linear(i, j, N)] for (i, j) in inds]
return np.all(np.equal(actual, expected))
"""Run `test` for a bunch of `N`s and make sure they're all right"""
print(all(map(test, range(2, 20))))
# prints True
值得写一篇博客 post 来解释如何得出这个结论,但现在就这样了。
三角矩阵泛化到任何维度。例如,假设 A[a,b,c,d,e] 仅当 0<=a<=b<=c<=d<=e< N 时才为非零。我们可以将 A 压缩为线性数组X 其中 A[a,b,c,d,e] = X[a+B[b]+C[c]+D[d]+E[e]]
其中B[b] = {b+1选2},C[c] = {c+2选3},D[d] = {d+3选4},E[e] = {e+4 选择 5}.
这些“偏移”数组可以在不使用乘法或除法的情况下计算如下:
B[0] = C[0] = D[0] = E[0] = 0;
for(int t = 1; t < N; t++)
{
B[t] = B[t-1]+t;
C[t] = C[t-1]+B[t];
D[t] = D[t-1]+C[t];
E[t] = E[t-1]+D[t];
}
我有一个矩阵的上三角部分,主对角线存储为线性数组,如何从数组的线性索引中提取矩阵元素的 (i,j) 索引?
例如线性数组:[a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10]
是矩阵的存储
a0 a1 a2 a3
0 a4 a5 a6
0 0 a7 a8
0 0 0 a10
我找到了这个问题的解决方案,但没有主对角线:
index = (n*(n-1)/2) - (n-i)*((n-i)-1)/2 + j - i - 1
同样问题的解决方案,但对于具有对角线的下三角矩阵:
index = ((i + 1) * i / 2 + i).
此致,
我找到答案了!它是:
index = (n*(n+1)/2) - (n-i)*((n-i)+1)/2 + j - i
我的解法可能等同于
index = N * i - ((i - 1) * i) / 2 + (j - i)
这是一个完整的 Python 测试。我使用 Python 因为 Numpy 有 triu_indices
,它给出了上三角索引。
import numpy as np
def mksquare(N):
"""Make a square N by N matrix containing 0 .. N*N-1"""
return np.arange(N * N).reshape(N, N)
def mkinds(N):
"""Return all triu indexes for N by N matrix"""
return [(i,j) for i in range(N) for j in range(N) if i <= j]
def ij2linear(i, j, N):
"""Convert (i,j) 2D index to linear triu index for N by N array"""
return N * i - ((i - 1) * i) // 2 + (j - i)
def test(N):
"""Make sure my `mkinds` works for given N"""
arr = mksquare(N)
vec = arr[np.triu_indices(N)]
inds = mkinds(N)
expected = [arr[i, j] for (i, j) in inds]
actual = [vec[ij2linear(i, j, N)] for (i, j) in inds]
return np.all(np.equal(actual, expected))
"""Run `test` for a bunch of `N`s and make sure they're all right"""
print(all(map(test, range(2, 20))))
# prints True
值得写一篇博客 post 来解释如何得出这个结论,但现在就这样了。
三角矩阵泛化到任何维度。例如,假设 A[a,b,c,d,e] 仅当 0<=a<=b<=c<=d<=e< N 时才为非零。我们可以将 A 压缩为线性数组X 其中 A[a,b,c,d,e] = X[a+B[b]+C[c]+D[d]+E[e]]
其中B[b] = {b+1选2},C[c] = {c+2选3},D[d] = {d+3选4},E[e] = {e+4 选择 5}.
这些“偏移”数组可以在不使用乘法或除法的情况下计算如下:
B[0] = C[0] = D[0] = E[0] = 0;
for(int t = 1; t < N; t++)
{
B[t] = B[t-1]+t;
C[t] = C[t-1]+B[t];
D[t] = D[t-1]+C[t];
E[t] = E[t-1]+D[t];
}