@[email protected] 源码分析: 数学后端-矩阵

[email protected]

这篇博客我们来看看更加基础层次上的关于PALISADE是如何实现常规数学运算的部分

关于为何要构建这样一个大型的数学库 我是这么想的

首先 第三方的库在PALISADE这样大规模工程环境中必定要做许多的源码调试和改动才能适应其需要的运算环境 但是另外一方面 这些库虽然可能更可靠，但是要直接上手去改现成的，那无疑是要对其开发的模式和约定有深刻的了解，并且很多功能可能是冗余的，在这种情况下，不如去做一套自己的小规模数学后端

一些经典的实现可能是必要的引用 但是在另一方面 针对格加密计算过程 数学后端可以给出更有针对性的优化策略

任何这样的加密系统想必都需要对一些常规的数据类型进行再封装，其中最有普遍性的就是矩阵和向量 让我们来palisade中对矩阵和向量的封装以及运算的实现

矩阵的乘运算:

template <class Element>

Matrix<Element> Matrix<Element>::Mult(Matrix<Element> const& other) const {

  // NUM_THREADS = omp_get_max_threads();

  if (cols != other.rows) {

    PALISADE_THROW(math_error, "incompatible matrix multiplication");

  }

  Matrix<Element> result(allocZero, rows, other.cols);

  if (rows == 1) {

#pragma omp parallel for

    for (size_t col = 0; col < result.cols; ++col) {

      for (size_t i = 0; i < cols; ++i) {

        result.data[0][col] += data[0][i] * other.data[i][col];

      }

    }

  } else {

#pragma omp parallel for

    for (size_t row = 0; row < result.rows; ++row) {

      for (size_t i = 0; i < cols; ++i) {

        for (size_t col = 0; col < result.cols; ++col) {

          result.data[row][col] += data[row][i] * other.data[i][col];

        }

      }

    }

  }

  return result;

}

可以看到矩阵的乘运算中 似乎是加入了对于并行(parallel)的优化设计

并不是以特定的并行计算的核分配方案，而是把这个任务交给了计算平台本身

矩阵计算乘是非常难以优化的 优化目标是朝着 n ^ 2 去逼近

这里可能有一个不那么现实的issue(事实上我也不知道有没有人提过了)，就是去调用开源的可行性优化方案

通常只针对一个问题的优化算法的引入并不会造成很大的嵌入负担 但是对于palisade这种应用广泛的系统性极强的同态加密库可能能由基础运算上带来可见的提高

参:

哈佛、MIT学者联手，创下矩阵乘法运算最快纪录 - 知乎 (zhihu.com)

2010.05846.pdf (arxiv.org)

// YSP The signature of this method needs to be changed in the future

// Laplace's formula is used to find the determinant

// Complexity is O(d!), where d is the dimension

// The determinant of a matrix is expressed in terms of its minors

// recursive implementation

// There are O(d^3) decomposition algorithms that can be implemented to support

// larger dimensions. Examples include the LU decomposition, the QR

// decomposition or the Cholesky decomposition(for positive definite matrices).

template <class Element>

void Matrix<Element>::Determinant(Element* determinant) const {

  if (rows != cols)

    PALISADE_THROW(math_error, "Supported only for square matrix");

  // auto determinant = *allocZero();

  if (rows < 1) PALISADE_THROW(math_error, "Dimension should be at least one");

  if (rows == 1) {

    *determinant = data[0][0];

  } else if (rows == 2) {

    *determinant = data[0][0] * (data[1][1]) - data[1][0] * (data[0][1]);

  } else {

    size_t j1, j2;

    size_t n = rows;

    Matrix<Element> result(allocZero, rows - 1, cols - 1);

    // for each column in sub-matrix

    for (j1 = 0; j1 < n; j1++) {

      // build sub-matrix with minor elements excluded

      for (size_t i = 1; i < n; i++) {

        j2 = 0;  // start at first sum-matrix column position

        // loop to copy source matrix less one column

        for (size_t j = 0; j < n; j++) {

          if (j == j1) continue;  // don't copy the minor column element

          // copy source element into new sub-matrix i-1 because new sub-matrix

          // is one row (and column) smaller with excluded minors

          result.data[i - 1][j2] = data[i][j];

          j2++;  // move to next sub-matrix column position

        }

      }

      auto tempDeterminant(allocZero());

      result.Determinant(&tempDeterminant);

      if (j1 % 2 == 0)

        *determinant = *determinant + (data[0][j1]) * tempDeterminant;

      else

        *determinant = *determinant - (data[0][j1]) * tempDeterminant;

      // if (j1 % 2 == 0)

      //  determinant = determinant + (*data[0][j1]) *

      // result.Determinant(); else   determinant = determinant -

      // (*data[0][j1]) * result.Determinant();

    }

  }

  // return determinant;

  return;

}

对于求行列式的值 显然实现上给出了更加丰富的注释 让我们先来读一下这段注释

YSP 这个方法以后需要改动，以拉普拉斯公式用于求行列式 复杂度是 O(d!)，其中 d 是维度 矩阵的行列式用其次要表示 递归实现 可以实现 O(d^3) 分解算法来支持 更大的可用性。示例包括 LU 分解、QR 分解或 Cholesky 分解（对于正定矩阵）。

从这段话来看 1. 作者应该有类似的优化企图 不过目前为止还没有实现(可能是因为维护团队内部的意见?) 2. 求行列式这个步骤作为一个可以被提高的部分其实可以作为往后优化的一个重要的目标 3. 这种优化可能是建立在对矩阵特征的分析和基于其导出的特征上的特定优化策略

而对于同样担负着重要角色的求逆工作 是这样实现的

// The cofactor matrix is the matrix of determinants of the minors A_{ij}

// multiplied by -1^{i+j} The determinant subroutine is used

template <class Element>

Matrix<Element> Matrix<Element>::CofactorMatrix() const {

  if (rows != cols)

    PALISADE_THROW(not_available_error, "Supported only for square matrix");

  size_t ii, jj, iNew, jNew;

  size_t n = rows;

  Matrix<Element> result(allocZero, rows, cols);

  for (size_t j = 0; j < n; j++) {

    for (size_t i = 0; i < n; i++) {

      Matrix<Element> c(allocZero, rows - 1, cols - 1);

      iNew = 0;

      for (ii = 0; ii < n; ii++) {

        if (ii == i) continue;

        jNew = 0;

        for (jj = 0; jj < n; jj++) {

          if (jj == j) continue;

          c.data[iNew][jNew] = data[ii][jj];

          jNew++;

        }

        iNew++;

      }

      Element determinant(allocZero());

      c.Determinant(&determinant);

      Element negDeterminant = -determinant;

      if ((i + j) % 2 == 0)

        result.data[i][j] = determinant;

      else

        result.data[i][j] = negDeterminant;

    }

  }

  return result;

}

逆矩阵还是一种非常常见的运算，如果使用平时手工计算工作使用的方式，严重影响了性能。而这里的计算方法则是我们平时经常使用的求代数余子式的方式。要在这一块进行优化确实是有相当难度的，因为关于求逆运算大多数也是因势制宜的，但是如果是对于随机生成的矩阵是没有什么特别通用的，大刀阔斧的优化方式的。

大多数场景下 优化的程度其实仍然和矩阵特性本身息息相关 比如使用高斯-约旦法去优化BillBoard矩阵，对于符合要求正定型采取特定的剪枝。这样的话，如果想要在这样的层次上去优化，就意味着增加一些判别开销来应对可能不太常见的情况。

这种情况下，如果是遇到可以被优化的矩阵类型则会大幅提高速度，但反过来，就会拖慢计算开始的时点进而导致计算速度有了微小的下移。毕竟能符合优化要求的矩阵出现概率不大，而在频繁的矩阵基础运算中，这可能是很可怕的，在较为复杂的计算流程(指该操作被重复调用的次数甚多)甚至会造成秒分级的速度下降。所以在这里我认同这里的实现方式，不采用特定的策略优化，而是去接收这种可以被允许的复杂计算。