将张量的"lower diagonal"映射到矩阵,作为将矩阵的下三角部分提取为向量的推广

Mapping the "lower diagonal" of a tensor to a matrix, as a generalization of the extraction of the lower triangular part of a matrix into a vector

给定一个四阶张量(每个阶的维度为K),例如T(p,q,r,s),我们可以将所有张量元素1对1映射到维度为[=的矩阵中14=],例如 M(i,j) 其中前两个张量索引 p,q 和最后两个索引 r,s 以列主要方式组合:

i = p + K * q
j = r + K * s

利用给定张量的一些(反)对称性,例如 T(p,q,r,s) = -T(q,p,r,s) = -T(p,q,s,r) = T(q,p,s,r)T(p,q,r,s) = T(r,s,p,q),我们希望能够构建一个矩阵 H(m,n),它只包含独特的元素(即那些与先前定义的对称性无关的元素),例如 p>qr>s 进入矩阵 H(m,n),然后维度为 K(K-1)/2 x K(K-1)/2.

我们如何找到一种算法(或者更好:我们如何使用 C++ 库 Eigen)来完成这些索引转换?此外,我们是否可以根据 p,qr,s 代数地写下 mn,就像我们在想要提取严格的下三角矩阵的情况下所做的那样矩阵(无对角线)转换为向量?

作为参考,给定一个方阵 Eigen::MatrixXd M (K,K),这是一个使用 C++ Eigen 库将给定方阵 M 的严格下三角提取到大小为 m 的向量的算法:

    Eigen::VectorXd strictLowerTriangle(const Eigen::MatrixXd& M) {

    auto K = static_cast<size_t>(M.cols());  // the dimension of the matrix
    Eigen::VectorXd m = Eigen::VectorXd::Zero((K*(K-1)/2));  // strictly lower triangle has K(K-1)/2 parameters

    size_t vector_index = 0;
    for (size_t q = 0; q < K; q++) {  // "column major" ordering for, so we do p first, then q
        for (size_t p = q+1; p < K; p++) {
            m(vector_index) = M(p,q);
            vector_index++;
        }
    }

    return m;
}

我们能够将此算法扩展到请求的一般情况:

Eigen::MatrixXd strictLowerTriangle(const Eigen::Tensor<double, 4>& T) {

    auto K = static_cast<size_t> (dims[0]);

    Eigen::MatrixXd M (K*(K-1)/2, K*(K-1)/2);

    size_t row_index = 0;
    for (size_t j = 0; j < K; j++) {  // "column major" ordering for row_index<-i,j so we do j first, then i
        for (size_t i = j+1; i < K; i++) {  // in column major indices, columns are contiguous, so the first of two indices changes more rapidly
                                                  // require i > j

            size_t column_index = 0;
            for (size_t l = 0; l < K; l++) {  // "column major" ordering for column_index<-k,l so we do l first, then k
                for (size_t k = l+1; k < K; k++) {  // in column major indices, columns are contiguous, so the first of two indices changes more rapidly
                                                          // require l > k

                    M(row_index,column_index) = T(i,j,k,l);

                    column_index++;
                }
            }

            row_index++;
        }
    }

    return M;
}