标签: matrix-multiplication

我正在寻找一个SIMD库,专注于小型(4x4)矩阵操作的图形.那里有很多单精度的,但我需要支持单精度和双精度.

我看过英特尔的IPP MX库,但我更喜欢有源代码的东西.我对这些特定操作的SSE3 +实现非常感兴趣:

1. Mat4*Mat4
2. Mat4*Vec4
3. Mat4*Mat4数组
4. Mat4*Vec4阵列
5. Mat4倒置(很高兴)

编辑:请不要"过早优化"答案.任何使用小矩阵的人都知道GCC不会对这些以及手动优化的内在函数或ASM进行矢量化.在这种情况下,这很重要,或者我不会问.

下面是C++实现,它比较了Eigen和For Loop执行矩阵 - 矩阵产品所用的时间.For循环已经过优化,可最大限度地减少缓存未命中.for循环最初比Eigen快,但最终变得更慢(500乘500矩阵可达2倍).我还应该做些什么才能与Eigen竞争？阻止了更好的本征性能的原因？如果是这样,我应该如何为for循环添加阻塞？

#include<iostream>
#include<Eigen/Dense>
#include<ctime>

int main(int argc, char* argv[]) {
    srand(time(NULL));
    //  Input the size of the matrix from the user
    int N   =   atoi(argv[1]);

    int M   =   N*N;

    //  The matrices stored as row-wise vectors
    double a[M];
    double b[M];
    double c[M];

    //  Initializing Eigen Matrices
    Eigen::MatrixXd a_E =   Eigen::MatrixXd::Random(N,N);
    Eigen::MatrixXd b_E =   Eigen::MatrixXd::Random(N,N);
    Eigen::MatrixXd c_E(N,N);

    double CPS  =   CLOCKS_PER_SEC;

    clock_t start, end;


    //  Matrix vector product by Eigen
    start   =   clock();
    c_E     =   a_E*b_E;
    end     = …

我知道两个完整矩阵乘法的下界是Ω(n ^ 2).矩阵乘法

我一直试图用问题转换方法证明两个下三角矩阵乘法的下界.

我最初的想法是(1)变换下三角矩阵,(2)估计这种变换的时间复杂度.

T(lower_triangular_matrix_multiplication(n))+O(lower_triangular_matrix_transformation(n))>?(full_matrix_multiplication(n)) = ?(n^2)

现在,我只需要证明O(lower_triangular_matrix_transformation(n)),我需要使三角矩阵成为一个完整的矩阵,所以为了简单起见,我只是让这个三角矩阵乘以它自身的变化,比如转置.

原因是下三角矩阵的平方仍然是下三角矩阵,下三角矩阵乘以其转置变化是"全矩阵".

所以我只需要分析三角矩阵的复杂性乘以其转置变化.

任何人都可以表明我的想法是否"合理"？

我在Android中有一个OpenGL ES View,由一个矩阵控制进行翻译.我试图找出一种方法来获得动态滚动的提示,如谷歌地图应用程序或iPhone中所见.谢谢.

我开发了一种科学应用(模拟在细胞核中移动的染色体).将染色体分成小片段,使用4x4旋转矩阵围绕随机轴旋转.

问题是模拟执行数千亿次旋转,因此浮点舍入误差会呈指数级增长并逐渐增长,因此随着时间的推移,碎片会"漂浮"并与染色体的其余部分分离.

我在C++中使用双精度.软件暂时在CPU上运行,但将移植到CUDA,模拟最多可持续1个月.

我不知道我怎么能以某种方式重新规范化染色体,因为所有的片段都被链接在一起(你可以把它看成是一个双重链接列表),但我认为如果可能的话,这将是最好的想法.

你有什么建议吗 ？我觉得有点迷茫.

非常感谢你,

H.

编辑:添加了简化的示例代码.您可以假设所有矩阵数学都是经典实现.

// Rotate 1000000 times
for (int i = 0; i < 1000000; ++i)
{
    // Pick a random section start
    int istart = rand() % chromosome->length;

    // Pick the end 20 segments further (cyclic)
    int iend = (istart + 20) % chromosome->length;

    // Build rotation axis
    Vector4 axis = chromosome->segments[istart].position - chromosome->segments[iend].position;
    axis.normalize();

    // Build rotation matrix and translation vector
    Matrix4 rotm(axis, rand() / float(RAND_MAX));
    Vector4 oldpos = chromosome->segments[istart].position; …

对于一个班级,我的老师提出的问题是将矩阵乘以其转置的算法成本.使用标准的3循环矩阵乘法算法,效率为O(N ^ 3),我想知道是否有办法操纵或利用矩阵*矩阵转置来获得更快的算法.我知道当你将矩阵乘以它的转置时,你必须计算较少的矩阵,因为它是对称的,但我想不出如何操作一个可能小于O(n ^ 3)的算法.

我知道像Coppensmith和Straussen这样的算法是更快的通用矩阵乘法算法,但任何人都可以提供有关如何计算利用转置的任何提示或见解？

谢谢

我有一个密钥算法,其中大部分运行时用于计算密集矩阵产品:

A*A'*Y, where: A is an m-by-n matrix, 
               A' is its conjugate transpose,
               Y is an m-by-k matrix

Typical characteristics:
    - k is much smaller than both m or n (k is typically < 10)
    - m in the range [500, 2000]
    - n in the range [100, 1000]

基于这些维度,根据矩阵链乘法问题的教训,很明显,在运算数意义上将计算结构化为最优A*(A'*Y).我当前的实现就是这样做的,而只是强迫关联性到表达式的性能提升是显而易见的.

我的应用程序是用C++编写的,用于x86_64平台.我正在使用Eigen线性代数库,英特尔的数学核心库作为后端.Eigen能够使用IMKL的BLAS接口来执行乘法,并且从我的Sandy Bridge机器上移动到Eigen的原生SSE2实现到Intel优化的基于AVX的实现的提升也很重要.

然而,表达式A * (A.adjoint() * Y)(用Eigen的说法)被分解为两个通用的矩阵 - 矩阵乘积(调用xGEMMBLAS例程),在它们之间创建一个临时矩阵.我想知道,通过一次专门的实现来一次评估整个表达式,我可以得到一个比我现在的通用更快的实现.一些让我相信这一点的观察结果是:

• 使用上述典型尺寸,输入矩阵A通常不适合缓存.因此,用于计算三矩阵乘积的特定存储器访问模式将是关键.显然,避免为部分产品创建临时矩阵也是有利的.

• A 并且它的共轭转置显然具有非常相关的结构,可以利用它来改善整体表达的存储器访问模式.

是否有任何标准技术以缓存友好的方式实现这种表达式？我发现的矩阵乘法的大多数优化技术都是针对标准A*B情况而不是更大的表达式.我对这个问题的微优化方面很满意,例如转换成适当的SIMD指令集,但是我正在寻找任何可以用尽可能最友好的方式打破这个结构的引用.

编辑: …

我正在执行NxN稀疏(~1-2%)矩阵的几个矩阵乘法,我们称之为B,使用NxM密集矩阵,我们称之为A(其中M <N).N很大,M也是如此; 大约几千.我正在运行Matlab 2013a.

现在,通常,矩阵乘法和大多数其他矩阵运算在Matlab中隐式并行化,即它们自动使用多个线程.如果任何一个矩阵都是稀疏的,那么看起来并非如此(参见例如StackOverflow讨论 - 没有对预期问题的答案 - 以及这个基本上没有答案的MathWorks线程).这对我来说是一个相当不愉快的惊喜.

我们可以通过以下代码验证多线程对稀疏矩阵操作没有影响:

clc; clear all; 

N = 5000;         % set matrix sizes
M = 3000;       
A = randn(N,M);   % create dense random matrices
B = sprand(N,N,0.015); % create sparse random matrix
Bf = full(B);     %create a dense form of the otherwise sparse matrix B

for i=1:3 % test for 1, 2, and 4 threads
  m(i) = 2^(i-1);
  maxNumCompThreads(m(i)); % set the thread count available to Matlab
  tic                      % starts …

我正在研究的算法需要在几个地方计算一种矩阵三元产品.

该操作采用具有相同尺寸的三个方形矩阵,并产生3指数张量.标记操作数A,B以及结果C的 (i,j,k)第th个元素

X[i,j,k] = \sum_a A[i,a] B[a,j] C[k,a]

在numpy中,你可以用它来计算einsum('ia,aj,ka->ijk', A, B, C).

问题:

• 此操作是否具有标准名称？
• 我可以通过一次BLAS呼叫来计算吗？
• 是否有任何其他重度优化的数值C/Fortran库可以计算这种类型的表达式？

我在图中有一组MxM 对称矩阵变量,其值我想优化.

有没有办法强制对称条件？

我已经考虑过为损失函数添加一个术语来强制执行它,但这看起来很尴尬和迂回.我希望的是tf.matmul(A,B,symmA=True) ,只有A的三角形部分才会被使用和学习.或者类似的东西tf.upperTriangularToFull(A)会从三角形部分创建一个密集的矩阵.