将稀疏矩阵乘以数组而不进行复制

Question

我有一个一维数组（大小为M的向量），一个很大的数组，我绝对不想在内存中复制它。我也有一个窗口N的稀疏矩阵（任意大小，基本上所有元素（对角线和N伪对角线除外）均为零）。

我想将此稀疏矩阵乘以向量，而不必将向量复制到内存中。最简单，最有效的方法是什么？必须有一个整洁的解决方案，但是我不了解合适的文献，而且我没有受过足够的教育，无法弄清楚这一点。

N=1有一个解决方案（其中矩阵是：对角线上的a，而两个最接近的伪对角线上的b）。解决方案看起来像这样（例如，在python中）：

tmp2 = array[0]
i = 1
while (i < len(array) - 1):
  tmp1 = b * (array[i - 1] + array[i + 1]) + a * array[i]
  array[i - 1] = tmp2
  i += 1
  tmp2 = b * (array[i - 1] + array[i + 1]) + a * array[i]
  array[i - 1] = tmp1
  i += 1

但是我无法将其概括为任意N。

Notes：我绝对不想在内存中复制大小为M的向量。但是，可以使用大小为2N+1的临时数组，因为M >> N。我正在寻找一个实际的算法描述，而不是一个能完成任务的智能自定义库。

提前感谢！

Answer 1

考虑矩阵

[
  1 2 3 0 0 0
  2 1 2 4 0 0
  3 2 1 2 5 0
  0 7 2 1 2 6
  0 0 8 2 1 2
  0 0 0 9 2 1
]

和向量v [1,2,3,4,5,6]

对于每一行，在v的涉及系数以下：

[1,2,3]
[1,2,3,4]
[1,2,3,4,5]
[  2,3,4,5,6]
[    3,4,5,6]
[      4,5,6]

您已经注意到，您只需要跟踪v的窗口。

该窗口最初是[1,2,3,4,5]（对于i = 0、1、2）

然后您将每个i窗口右移（并最终将其截断以使最后一行不超出v的范围...]

现在注意，当您向右移动时，您只需要知道v中的下一个值，只要您没有弄脏该值（通过写入v），您的新窗口就有效。] >

对于行i，窗口为[i-n;i+n]，将被修改的系数为v[i]。对于下一个窗口，您需要知道尚未被弄脏的v[i+n+1]。一切都很好。

所以总是像

window = circularbuffer(2n+1) //you push to the right, and if length > 2n+1, you remove first elem
for i = 0; i<v.size()
  v[i] = prod(row_i, window) // only for the row_i coeffs...
  if i >= n && < M-3
    window.push(v[i+n+1])
  else if i>= M-3
    window.shift() // just remove the first value

const N = 2
const M_SIZE = 10
function toString(M){
  return M.map(x=>x.join(' ')).join('\n')
}
const { M, v } = (_ => {
  const M = Array(M_SIZE).fill(0).map(x=>Array(M_SIZE).fill(0))
  let z = 1
  for(let i = 0; i<M_SIZE; ++i){
    for(let j = -N; j<=N; ++j){
      if(i+j >= 0 && i+j <M_SIZE){
        M[i][i+j] = (z++ % (N*2))+1
      }
    }
  }
  const v = Array(M.length).fill(0).map((x,i)=>i)
  return { M, v}
})()

function classic(M, v){
  return M.map(r => r.reduce((acc, x, j) => acc + v[j]*x, 0))
}

function inplace(M, v){
  // captn inefficiency
  const circBuf = (init => {
    let buf = init
    return {
      push (x) {
        buf.push(x)
        buf.shift()
      },
      shift() {
        buf.shift()
      },
      at (i) { return buf[i] },
      toString() {
        return buf.join(' ')
      }
    }
  })(v.slice(0, 2 * N + 1))

  const sparseProd = (row, buf) => {
    let s = 0
    row.forEach((x, j) => s += x * buf.at(j))
    return s
  }

  const sparseRows = M.map(r => r.filter(x => x !== 0))

  sparseRows.forEach((row, i) => {
    v[i] = sparseProd(row, circBuf)
    if (i >= sparseRows.length - 3 ) {
      circBuf.shift()
    } else {
      if (i >= N) {
        circBuf.push(v[i + N + 1])
      } 
    }
  })
}
console.log('classic prod', classic(M, v))
inplace(M, v)
console.log('inplace prod', v)