三角矩阵系数的索引号算法答案

【问题标题】：algorithm for index numbers of triangular matrix coefficients三角矩阵系数的索引号算法
【发布时间】：2008-10-28 09:58:14
【问题描述】：

我认为这一定很简单，但我无法正确...

我有一个 MxM 三角矩阵，其系数逐行存储在一个向量中。例如：

M =   [ m00 m01 m02 m03 ] 
      [     m11 m12 m13 ]
      [         m22 m23 ]
      [             m33 ]

存储为

coef[ m00 m01 m02 m03 m11 m12 m13 m22 m23 m33 ]

现在我正在寻找一种非递归算法，它可以为我提供矩阵大小M 和系数数组索引i

unsigned int row_index(i,M)

和

unsigned int column_index(i,M)

它引用的矩阵元素。所以， row_index(9,4) == 3、column_index(7,4) == 2 等，如果索引计数是从零开始的。

编辑：已经给出了几个使用迭代的回复。有人知道代数表达式吗？

【问题讨论】：

您能否改写一下，以便数组中的元素可能是字母而不是数字 - 当您同时使用零时，很难 100% 准确理解您的函数应该做什么数组本身中基于行/列的值和基于 1 的值。
data[row * num_cols + column] 是作为单个数组存储在内存中的 C 二维数组的经典表达式。
显然矩阵中重复的m12是错误的。
是的，第二个应该是 m13。我想在下面的解决方案中提到它，但没有:-)

标签： algorithm

【解决方案1】：

在此回复的最后一行，解释如下:-)

系数数组具有：第一行（第 0 行，在您的索引中）的 M 个元素，第二行（第 1 行）的 (M-1) 个元素，依此类推，总共 M+(M-1) +…+1=M(M+1)/2 个元素。

从最后开始工作稍微容易一些，因为系数数组总是在最后一行（第 M-1 行）有 1 个元素，在倒数第二行（第 M- 行）有 2 个元素2)，第 M-3 行的 3 个元素，依此类推。最后K行占据了系数数组的最后K(K+1)/2个位置。

现在假设给定系数数组中的索引 i。在大于 i 的位置有 M(M+1)/2-1-i 个元素。拨打这个号码 i'; 你想找到最大的 k 使得 k(k+1)/2 ≤ i'。这个问题的形式是你痛苦的根源——据我所知，你无法避免取平方根:-)

不管怎样，让我们这样做：k(k+1) ≤ 2i' 意味着 (k+1/2)² - 1/4 ≤ 2i'，或者等效地 k ≤ (sqrt(8i '+1)-1)/2。让我把最大的 k 称为 K，那么后面有 K 行（总共 M 行），所以 row_index(i,M) 是 M-1-K。至于列索引，i'中的K(K+1)/2个元素在后面的行中，所以这一行后面的元素有j'=i'-K(K+1)/2个（有总共K+1个元素），所以列索引是K-j'。 [或者等效地，这一行从末尾开始是(K+1)(K+2)/2，所以我们只需要从i中减去这一行的起始位置：i-[M(M+1)/2 -(K+1)(K+2)/2]。令人振奋的是，这两种表达方式都给出了相同的答案！]

（伪）代码[使用 ii 而不是 i'，因为某些语言不允许名称像 i' ;-)]：

row_index(i, M):
    ii = M(M+1)/2-1-i
    K = floor((sqrt(8ii+1)-1)/2)
    return M-1-K

column_index(i, M):
    ii = M(M+1)/2-1-i
    K = floor((sqrt(8ii+1)-1)/2)
    return i - M(M+1)/2 + (K+1)(K+2)/2

当然，您可以通过将表达式替换 ii 和 K 将它们变成单行。您可能需要注意精度误差，但有一些方法可以找到整数平方根，而无需浮点计算.另外，引用 Knuth 的话：“当心上述代码中的错误；我只是证明它是正确的，没有尝试过。”

如果我可以冒昧地进一步说明一下：简单地将所有值保存在 M×M 数组中最多会占用两倍的内存，并且根据您的问题，与算法改进相比，2 的因子可能微不足道，并且可能值得用可能昂贵的平方根计算来换取更简单的表达式。

[编辑：顺便说一句，您可以证明 floor((sqrt(8ii+1)-1)/2) = (isqrt(8ii+1)-1)/2 其中 isqrt(x)=floor(sqrt( x)) 是整数平方根，除法是整数除法（截断；C/C++/Java 等中的默认值）——所以如果你担心精度问题，你只需要担心实现一个正确的整数平方根。]

【讨论】：

这正是我需要找到的！就我而言，我试图对 GPU 上两个向量之间的每一对元素执行计算，因此行和列是元素索引。在这种情况下，它不仅仅是将内存占用量增加一倍，尽管这也是问题的很大一部分。我需要做一些调整来处理内存寻址，但你的解释给出了一个很好的框架。谢谢！
你能展示一下没有对角元素的三角矩阵（即 M(M-1)/2 个元素）的代数是什么样子的吗？
@frank 当然，它是相似的。最后 K 行共有 K(K-1)/2 个元素。所以给定 i（基于 0 的索引），我们想要最大的 K 使得 K(K-1)/2 ≤ M(M-1)/2 - 1 - i。这适用于 K = floor(1/2 * (1 + sqrt(4M*(M-1) - 8*i + 7)))。使用该 K 值，行索引（从 0 开始）为 M-1-K。
我不确定答案是否有小错误。此答案的列索引似乎从对角元素开始，而不是完整矩阵中的列索引。你能验证一下吗？

【解决方案2】：

这是一个代数（主要）解决方案：

unsigned int row_index( unsigned int i, unsigned int M ){
    double m = M;
    double row = (-2*m - 1 + sqrt( (4*m*(m+1) - 8*(double)i - 7) )) / -2;
    if( row == (double)(int) row ) row -= 1;
    return (unsigned int) row;
}


unsigned int column_index( unsigned int i, unsigned int M ){
    unsigned int row = row_index( i, M);
    return  i - M * row + row*(row+1) / 2;
}

编辑：固定 row_index()

【讨论】：

我不认为这适用于每个 M。例如，该解决方案给出 row_index(12,5)=2，而正确答案（如果我做得正确的话）是 3。
我的观察也是如此。 M=6 为 row_index 生成 3 个错误值（因此也为 column_index 生成）。

【解决方案3】：

ShreevatsaR 的解释非常好，帮助我解决了我的问题。但是，为列索引提供的解释和代码给出的答案与问题要求的答案略有不同。

重申一下，在 i 之后的行中有 j' = i' - K(K+1)/2 个元素。但是该行与其他行一样，有 M 个元素。因此（从零开始的）列索引是 y = M - 1 - j'。

对应的伪代码为：

column_index(i, M):
  ii = M(M+1)/2-1-i
  K = floor((sqrt(8ii+1)-1)/2)
  jj = ii - K(K+1)/2
  return M-1-jj

ShreevatsaR 给出的答案，K - j'，是从矩阵的对角线开始计数（为零）时的列索引。因此，他的计算得出 column_index(7,4) = 0 而不是问题中指定的 column_index(7,4) = 2。

【讨论】：

你能展示一下没有对角元素的三角矩阵的代数是什么样的吗（即对于 M(M-1)/2 个元素，索引也从 0 开始）？

【解决方案4】：

这些可能有一个聪明的单行，但是（减去任何错误检查）：

unsigned int row_index( unsigned int i, unsigned int M ){
    unsigned int row = 0;
    unsigned int delta = M - 1;
    for( unsigned int x = delta; x < i; x += delta-- ){
        row++;
    }
    return row;
}

unsigned int column_index( unsigned int i, unsigned int M ){
    unsigned int row = 0;
    unsigned int delta = M - 1;
    unsigned int x;
    for( x = delta; x < i; x += delta-- ){
        row++;
    }
    return M + i - x - 1;
}

【讨论】：

行索引OK，列索引偏移量为2，能不能改一下return M + i - x - i;的return语句？谢谢。

【解决方案5】：

应该是这样的

i == col + row*(M-1)-row*(row-1)/2

所以找到 col 和 row 的一种方法是遍历可能的 row 值：

for(row = 0; row < M; row++){
  col = i - row*(M-1)-row*(row-1)/2
  if (row <= col < M) return (row,column);
}

这至少是非递归的，我不知道你是否可以不迭代。

从这个答案和其他答案可以看出，您几乎必须计算行才能计算列，因此在一个函数中同时完成这两项工作可能是明智的。

【讨论】：

【解决方案6】：

我想了一下，得到了以下结果。请注意，您可以在一次拍摄中同时获得行和列。

假设：行从 0 开始。列从 0 开始。索引从 0 开始

符号

N = 矩阵大小（在原始问题中是 M）

m = 元素的索引

伪代码是

function ind2subTriu(m,N)
{
  d = 0;
  i = -1;
  while d < m
  {
    i = i + 1
    d = i*(N-1) - i*(i-1)/2
  }
  i0 = i-1;
  j0 = m - i0*(N-1) + i0*(i0-1)/2 + i0 + 1;
  return i0,j0
}

还有一些 octave/matlab 代码

function [i0 j0]= ind2subTriu(m,N)
 I = 0:N-2;
 d = I*(N-1)-I.*(I-1)/2;
 i0 = I(find (d < m,1,'last'));
 j0 = m - d(i0+1) + i0 + 1;

你怎么看？

截至 2011 年 12 月，执行此操作的非常好的代码已添加到 GNU/Octave。他们可能会扩展 sub2ind 和 ind2sub。该代码暂时可以作为私有函数 ind2sub_tril 和 sub2ind_tril 找到

【讨论】：

娟皮：你好像有一个旧的，未注册的帐户，和一个新的注册帐户。我已经标记了帖子以引起模组的注意，但您应该在 meta 上发帖要求合并帐户。
@JuanPi - 标记此问题并注明您的两个帐户，我们很乐意将其理顺。

【解决方案7】：

我花了一些时间来了解您需要什么！ :)

unsigned int row_index(int i, int m)
{
    int iCurrentRow = 0;
    int iTotalItems = 0;
    for(int j = m; j > 0; j--)
    {
        iTotalItems += j;

        if( (i+1) <= iTotalItems)
            return iCurrentRow;

        iCurrentRow ++;
    }

    return -1; // Not checking if "i" can be in a MxM matrix.
}

抱歉忘记了其他功能.....

unsigned int column_index(int i, int m)
{
    int iTotalItems = 0;
    for(int j = m; j > 0; j--)
    {
        iTotalItems += j;

        if( (i+1) <= iTotalItems)
            return m - (iTotalItems - i);
    }

    return -1; // Not checking if "i" can be in a MxM matrix.
}

【讨论】：