标签: matrix-multiplication

我想知道是否有人可以告诉我如何有效地使用循环平铺/循环阻塞来进行大密集矩阵乘法.我正在使用1000x1000矩阵进行C = AB.我已经按照维基百科上的示例进行循环平铺,但是使用平铺比使用平铺更糟糕.

http://en.wikipedia.org/wiki/Loop_tiling

http://software.intel.com/en-us/articles/how-to-use-loop-blocking-to-optimize-memory-use-on-32-bit-intel-architecture

我在下面提供了一些代码.由于缓存未命中,天真的方法非常慢.转置方法在缓冲区中创建B的转置.这种方法给出了最快的结果(矩阵乘法为O(n ^ 3)并且转置为O(n ^ 2),因此转置至少快1000倍).没有阻塞的wiki方法也很快,不需要缓冲区.阻塞方法较慢.阻塞的另一个问题是它必须多次更新块.这对线程/ OpenMP来说是一个挑战,因为如果不小心它会导致竞争条件.

我应该指出,使用AVX修改转置方法,我得到的结果比Eigen快.但是,我的SSE结果比Eigen慢一点,所以我想我可以更好地使用缓存.

编辑:我想我知道自己想做什么.它来自1969年的Cannon算法.http:
//en.wikipedia.org/wiki/Matrix_multiplication#Communication-avoiding_and_distributed_algorithms

我需要做的块矩阵乘法,并具有每个线程处理的一个子矩阵C ^而非阿和乙.例如,如果我将矩阵分成四个块.我会做:

+-+      +-+     +-+      +-+   +-+      +-+
|          |     |          |   |          |
| C1    C2 |     | A1    A2 |   | B1    B2 |
|          |  =  |          | x |          |
| C3    C4 |     | A3    A4 |   | B3    B4 |
|          |     |          |   |          |
+-+      +-+     +-+      +-+ …

这里使用hdf5的矩阵乘法我使用hdf5(pytables)进行大矩阵乘法,但我很惊讶因为使用hdf5它的工作速度更快,然后在RAM中使用普通的numpy.dot和存储矩阵,这种行为的原因是什么？

也许在python中有一些更快的矩阵乘法函数,因为我仍然使用numpy.dot进行小块矩阵乘法.

这是一些代码:

假设矩阵可以适合RAM:在矩阵10*1000 x 1000上进行测试.

使用默认numpy(我认为没有BLAS lib).普通的numpy数组在RAM中:时间9.48

如果A,B在RAM中,C在磁盘上:时间1.48

如果磁盘上的A,B,C:时间372.25

如果我使用numpy与MKL结果是:0.15,0.45,43.5.

结果看起来很合理,但我仍然不明白为什么在第一种情况下块乘法更快(当我们将A,B存储在RAM中时).

n_row=1000
n_col=1000
n_batch=10

def test_plain_numpy():
    A=np.random.rand(n_row,n_col)# float by default?
    B=np.random.rand(n_col,n_row)
    t0= time.time()
    res= np.dot(A,B)
    print (time.time()-t0)

#A,B in RAM, C on disk
def test_hdf5_ram():
    rows = n_row
    cols = n_col
    batches = n_batch

    #using numpy array
    A=np.random.rand(n_row,n_col)
    B=np.random.rand(n_col,n_row)

    #settings for all hdf5 files
    atom = tables.Float32Atom() #if store uint8 less memory?
    filters = tables.Filters(complevel=9, complib='blosc') # tune parameters
    Nchunk = 128  # ?
    chunkshape = (Nchunk, …

我不太明白是什么使C#/ .NET(甚至Java)中的矩阵乘法变得如此之慢.

看一下这个基准测试(来源):试图找到更新的基准测试.

Java vs C#vs C++细分http://img411.imageshack.us/img411/9324/perf.gif

C#的整数和双重性能非常接近用MSVC++编译的C++.双精度为87%,32位整数为99%.非常好,我会说.但接下来看矩阵乘法.差距扩大到C#的速度约为19%.这是一个非常大的差异,我不明白.矩阵乘法只是一堆简单的数学运算.怎么这么慢？它不应该与等量的简单浮点或整数运算一样快吗？

这尤其是游戏和XNA的关注点,其中矩阵和矢量性能对物理引擎等事物至关重要.前一段时间,Mono通过一些漂亮的矢量和矩阵类增加了对SIMD指令的支持.它缩小了差距,使Mono比手写的C++更快,尽管没有C++与SIMD一样快.(来源)

矩阵乘法比较http://img237.imageshack.us/img237/2788/resultse.png

这里发生了什么？

编辑:仔细观察,我误读了第二张图.C#看起来非常接近.第一个基准测试只是做了可怕的错误吗？对不起,我错过了第一个基准测试的版本号.我抓住它作为我总是听到的"C#线性代数很慢"的方便参考.我会试着找另一个.

在Haskell中写这个Scala矩阵乘法的冗长讨论后,我想知道......类型安全矩阵乘法会是什么样的？所以这是你的挑战:要么链接到Haskell实现,要么自己实现,如下:

data Matrix ... = ...

matrixMult :: Matrix ... -> Matrix ... -> Matrix ...
matrixMult ... = ...

当matrixMult产生一个错误类型在编译的时候,如果你试图乘两个matricies不兼容的尺寸.布朗尼指出,如果你链接到讨论这个精确主题的论文或书籍,和/或讨论自己这个功能是多么有用/无用.

我试图使用SSE找到使用向量(u)的4x4矩阵(M)乘法的最有效实现.我的意思是Mu = v.

据我所知,有两种主要方法可以解决这个问题:

    method 1) v1 = dot(row1, u), v2 = dot(row2, u), v3 = dot(row3, u), v4 = dot(row4, u)
    method 2) v = u1 col1 + u2 col2 + u3 col3 + u4 col4.

方法2易于在SSE2中实现.方法1可以用SSE3中的水平加法指令或SSE4中的点积指令来实现.但是,在我的所有测试中,方法2总是优于方法1.

我认为方法1有优势的一个地方是3x4矩阵,例如仿射变换.在这种情况下,最后一个点积是不必要的.但即使在这种情况下,4x4矩阵上的方法2也比3x4矩阵上的方法1快.我发现的唯一方法比4x4矩阵上的方法2快4x3矩阵上的方法2.

那么水平加法和点积指令有什么意义呢？事实上,点生产指令在这种情况下表现最差.也许它与数据格式有关？如果无法定义矩阵的排序方式,那么转置是必要的,在这种情况下,方法1可能会更好吗？

请参阅下面的一些代码.

__m128 m4x4v_colSSE(const __m128 cols[4], const __m128 v) {
  __m128 u1 = _mm_shuffle_ps(v,v, _MM_SHUFFLE(0,0,0,0));
  __m128 u2 = _mm_shuffle_ps(v,v, _MM_SHUFFLE(1,1,1,1));
  __m128 u3 = _mm_shuffle_ps(v,v, _MM_SHUFFLE(2,2,2,2));
  __m128 u4 = _mm_shuffle_ps(v,v, _MM_SHUFFLE(3,3,3,3));

  __m128 prod1 = _mm_mul_ps(u1, cols[0]);
  __m128 prod2 = _mm_mul_ps(u2, cols[1]);
  __m128 …

我正在寻找一种更快,更棘手的方法来将C中的两个4x4矩阵相乘.我目前的研究主要集中在具有SIMD扩展的x86-64汇编上.到目前为止,我已经创建了一个函数,比一个简单的C实现快6倍,这超出了我对性能改进的期望.不幸的是,只有在没有使用优化标志进行编译时(GCC 4.7),这种情况才会成立.随着-O2,C变得更快,我的努力变得毫无意义.

我知道现代编译器利用复杂的优化技术来实现几乎完美的代码,通常比巧妙的手工装配更快.但在少数性能关键的情况下,人类可能会尝试使用编译器争取时钟周期.特别是,当一些支持现代ISA的数学可以被探索时(就像我的情况一样).

我的函数如下(AT&T语法,GNU汇编程序):

    .text
    .globl matrixMultiplyASM
    .type matrixMultiplyASM, @function
matrixMultiplyASM:
    movaps   (%rdi), %xmm0    # fetch the first matrix (use four registers)
    movaps 16(%rdi), %xmm1
    movaps 32(%rdi), %xmm2
    movaps 48(%rdi), %xmm3
    xorq %rcx, %rcx           # reset (forward) loop iterator
.ROW:
    movss (%rsi), %xmm4       # Compute four values (one row) in parallel:
    shufps $0x0, %xmm4, %xmm4 # 4x 4FP mul's, 3x 4FP add's 6x mov's per row,
    mulps %xmm0, %xmm4        # expressed in four sequences of 5 instructions,
    movaps …

Octave有没有办法只计算和存储矩阵乘积的对角线？

基本上喜欢这样做: vector = diag(A*B);

我不关心A*B除了对角线上的任何值.矩阵大小是围绕80k x 12和12 x 80k,所以即使我不关心速度/额外的内存它只是惯于适合在RAM中.

奇怪的是,因为Octave是一个包含大量数据集和对角线的包非常重要,所以它应该是可能的.

我只需要矩阵乘法中的对角线元素:

,

在R.由于Z很大,我想避免全面的乘法......

Z <- matrix(c(1,1,1,2,3,4), ncol = 2)
Z
#     [,1] [,2]
#[1,]    1    2
#[2,]    1    3
#[3,]    1    4

X <- matrix(c(10,-5,-5,20), ncol = 2)
X
#     [,1] [,2]
#[1,]   10   -5
#[2,]   -5   20

Z %*% D %*% t(Z)
#     [,1] [,2] [,3]
#[1,]   70  105  140
#[2,]  105  160  215
#[3,]  140  215  290

diag(Z %*% D %*% t(Z))
#[1]  70 160 290

X始终是一个小方阵(2x2,3x3或4x4),其中Z的列数等于X的维数.是否有可用的函数？

假设我们有一个单一频道图像(5x5)

A = [ 1 2 3 4 5
      6 7 8 9 2
      1 4 5 6 3
      4 5 6 7 4
      3 4 5 6 2 ]

和过滤器K(2x2)

K = [ 1 1
      1 1 ]

应用卷积的一个例子(让我们从A中取出第一个2x2)将是

1*1 + 2*1 + 6*1 + 7*1 = 16

这非常简单.但是,让我们向矩阵A引入深度因子,即在深度网络中具有3个通道或甚至转换层的RGB图像(深度= 512).如何使用相同的滤波器完成卷积运算？类似的工作对于RGB情况非常有帮助.

我有2矩阵100kx200和200x100k如果它们是小矩阵我会使用numpy dot产品

sum(a.dot(b), axis = 0)

但矩阵太大了,我也不能使用循环是否有一个聪明的方法来做到这一点？