Отображение «нижней диагонали» тензора на матрицу как обобщение извлечения нижней треугольной части матрицы в вектор - PullRequest
Новая
       7

Отображение «нижней диагонали» тензора на матрицу как обобщение извлечения нижней треугольной части матрицы в вектор

0 голосов
/ 07 июня 2018

Учитывая тензор ранга 4 (каждый ранг с размерностью K), например, T(p,q,r,s), мы можем отобразить 1-к-1 все элементы тензора в матрицу размерности K^2 x K^2, например, * 1004.* в котором два первых индекса тензора 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>q и r>s вматрица H(m,n), которая тогда будет иметь размерность K(K-1)/2 x K(K-1)/2.

Как мы можем найти алгоритм (или даже лучше: как мы можем использовать библиотеку C ++ Eigen) для выполнения этих преобразований индекса?Кроме того, можем ли мы записать m и n алгебраически в терминах p,q и r,s, как мы можем это сделать в случае, когда мы хотим извлечь строгую нижнюю треугольную матрицу (без диагонали) в вектор

1 Ответ

0 голосов
/ 07 июня 2018

Для справки: для квадратной матрицы Eigen::MatrixXd M (K,K) приведен алгоритм, который выделяет строгий нижний треугольник данной квадратной матрицы M в вектор размером m с использованием библиотеки C ++ Eigen: 

    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;
}
...