java中的矩阵运算-matlab 矩阵或运算

选自TowardsDataScience

作者：Anuradha Wickramarachchi

机器之心编译

参与：Nurhachu Null

在机器学习中，绝大多数任务会涉及到耗费时间的大量运算，而且随着数据集的增加，运算量会越来越大。解决这个问题的一个方法就是使用多线程。在这篇文章中，我要结合代码介绍一下 GPU 加速，它是如何完成的，以及用于 GPU 任务的简单 API。下面以一个矩阵乘法开始全文内容。

矩阵乘法

上面给出了两个矩阵java中的矩阵运算，一个 3×6 的，一个 6×6 的。乘积的结果将会是一个 3×6 的矩阵。完成这个运算总共需要 3×6×6 次乘法运算。那么，我们可以得到这样的结论：这个任务的时间复杂度是 O（mn^2）。这也就意味着，2000×2000 的矩阵运算将会需要 8,000,000,000 次乘法运算。这会花费大量的 CPU 计算时间。

引入 GPU

通常 GPU 会包含大量的处理核心。核心数目从 384 个到几千个。下面是 NVIDIA 几款消费级 GPU 的比较（）：

CUDA 核数目

CUDA 是统一计算设备架构（Compute Unified Device Architecture）的缩写。它们以相对稍慢的速度运行，但是能够通过使用大量运算逻辑单元（ALU）来提供很大的并行度。更详细内容请参考链接（）。

CUDA 线程模型

这张图展示了 CUDA 的线程模型（这个和市场上其他的架构几乎是相同的，例如 AMD）。简单起见，我们假设一每个 CUDA 核一次只能运行一个线程。如果我们的数据集比较大，我们可以将它分成块。上图中的一个 Grid 包含多个 Block。Block 则是另一个包含与它维度相同个数的线程的矩阵。总之，由于这是一个简介，所以我们要以一个用 Java 开发的简单 API 来聚焦更大更复杂的结构。

GPU 的思考

正如我们讨论到的，每个 GPU 核心都能运行一个独立的线程。开始这个模拟的最简单的方式就是假设最终结果数组中的每个元素都由一个 GPU 核来计算。因为所有的核都是并行运行的，所有矩阵的所有元素也会被并行的计算。所以，我们现在的时间复杂度就变成了 O(n)。现在，对于 2000×2000 的矩阵乘法，我们只需要 2000 次运行，这对计算机而言是容易计算的。通常我们之前所说的每一个线程都知道自己的身份，也就是它所属于的 block 和 Grid。或者，说得简单一些就是元素在矩阵中的位置。此外，矩阵会被加载到 GPU 中共享它的内存，我们可以通过索引直接访问元组中的数据。是不是很容易？我们对着代码来看一看吧。

使用 APARAPI 进行 GPU 编程

APARAPI（A-PARallel-API）是一个基于 OpenCL 的用于 GPU 编程的 wrapper。它既支持 CUDA 架构，也支持 AMD 架构。此外，这个 API 还引入了 Java 中的伟大的面向对象思想，如果我们直接用 C++来完成这个任务的话也许会有些混乱。上手非常容易。虽然其中有内在依赖项，但是要确保你正确地设置了 OpenCL 或者 CUDA。简单的 Google 一下会帮助到你。大多数设备都是自带的（OSX 或者 windows 设备）。

pom.xml

xml version="1.0" encoding="UTF-8"?>
<project xmlns="http://maven.apache.org/POM/4.0.0"
         xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
         xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/xsd/maven-4.0.0.xsd">
    <modelVersion>4.0.0modelVersion>
    <groupId>cuda-aparapigroupId>
    <artifactId>cudaartifactId>
    <version>1.0-SNAPSHOTversion>
    <dependencies>
        <dependency>
            <groupId>com.aparapigroupId>
            <artifactId>aparapiartifactId>
            <version>1.4.1version>
        dependency>
    dependencies>
project>

MatrixMultiplication.java

import com.aparapi.Kernel;
import com.aparapi.Range;
/**
 * Created by anuradhawick on 12/29/17.
 */
public class MatrixMultiplication {
    public static void main(String[] args) {
        // Width of the matrix
        final int SIZE = 5000;
        // We should use linear arrays as supported by the API
        final int[] a = new int[SIZE * SIZE];
        final int[] b = new int[SIZE * SIZE];
        int[] c = new int[SIZE * SIZE];
        final int[] d = new int[SIZE * SIZE];
        int val;
        // Creating random matrices
        for (int i = 0; i < SIZE; i++) {
            for (int j = 0; j < SIZE; j++) {
                a[i * SIZE + j] = (int) (Math.random() * 100);
                b[i * SIZE + j] = (int) (Math.random() * 100);
            }
        }
        long time = System.currentTimeMillis();
        // CPU multiplication
        System.out.println("Starting single threaded computation");
        for (int i = 0; i < SIZE; i++) {
            for (int j = 0; j < SIZE; j++) {
                val = 0;
                for (int k = 0; k < SIZE; k++) {
                    val += a[i * SIZE + k] * b[k * SIZE + j];
                }
                c[i * SIZE + j] = val;
            }
        }
        System.out.println("Task finished in " + (System.currentTimeMillis() - time) + "ms");
        // Kernel for multiplication
        Kernel kernel = new Kernel() {
            @Override
            public void run() {
                int row = getGlobalId() / SIZE;
                int col = getGlobalId() % SIZE;
                if (row > SIZE || col > SIZE) return;
                d[row * SIZE + col] = 0;
                for (int i = 0; i < SIZE; i++) {
                    d[row * SIZE + col] += a[row * SIZE + i] * b[i * SIZE + col];
                }
            }
        };
        // Array size for GPU to know
        Range range = Range.create(SIZE * SIZE);
        System.out.println("Starting GPU computation");
        time = System.currentTimeMillis();
        kernel.execute(range); // Running the Kernel
        System.out.println("Task finished in " + (System.currentTimeMillis() - time) + "ms");
        // Verifying the result
        for (int i = 0; i < SIZE; i++) {
            for (int j = 0; j < SIZE; j++) {
                if (c[i * SIZE + j] != d[i * SIZE + j]) {
                    System.out.println("ERROR");
                    return;
                }
            }
        }
    }
}

上述代码的精简化

Kernel 就是在 GPU 上运行的代码部分。Kernel 可见的变量将会被拷贝到 GPU 的 RAM 中。我们因为 GPU 支持线性数组，所以我们不能以 2D 数组的形式输入数据。GPU 不能处理 2D 数组，但是它们是通过维度的概念来处理的（此处暂且不讨论这个内容）。

Range range = Range.create(SIZE * SIZE);

上述代码在 GPU 中分配了小于等于 SIZE × SIZE 个线程。

int row = getGlobalId() / SIZE;
int col = getGlobalId() % SIZE;

上述代码从私有内存中得到了线程的 ID。我们可以通过这个 ID 来区分这个线程单元的位置。对每个线程我们做以下处理：

for (int i = 0; i < SIZE; i++) {
    d[row * SIZE + col] += a[row * SIZE + i] * b[i * SIZE + col];
}

这是两个矩阵对应单元相乘相加的最简单的形式。我们只为使用线程索引的单个线程定义了 Kernel，它将会在所有的线程上并行运行。

结果

运算是很快的，但是有多快呢？这个是上述代码的输出：

1200 × 1200

Starting single threaded computation
Task finished in 25269ms 
Starting GPU computation
Task finished in 1535ms

由于下面的矩阵比较大java中的矩阵运算，所以我们只在 GPU 上运行以下的运算。

2000 × 2000 的矩阵运算耗时 3757ms

5000 × 5000 的矩阵运算耗时 5402ms

自己也尝试一下？

原文链接：