优于立方复杂度的 Rust 中矩阵乘法

中途:三次矩阵乘法

一、说明

        几年前,我在 C++ 年编写了 Strassen 矩阵乘法算法的实现,最近在 Rust 中重新实现了它,因为我继续学习该语言。这是学习 Rust 性能特征和优化技术的有用练习,因为尽管 Strassen 的算法复杂性优于朴素方法,但它在算法结构中的分配和递归开销中具有很高的常数因子

  • 通用算法
  • 换位以获得更好的性能
  • 次立方:斯特拉森算法的工作原理
  • 排比
  • 标杆
  • 分析和性能优化

二、通用算法

        一般(朴素)矩阵乘法算法是每个人在他们的第一堂线性代数课上学习的三个嵌套循环方法,大多数人会将其识别为 O(n³)

pub fn 
mult_naive (a: &Matrix, b: &Matrix) -> Matrix {if a.rows == b.cols {let m = a.rows;let n = a.cols;// preallocatelet mut c: Vec<f64> = Vec::with_capacity(m * m);for i in 0..m {for j in 0..m {let mut sum: f64 = 0.0;for k in 0..n {sum += a.at(i, k) * b.at(k, j);}c.push(sum);}}return Matrix::with_vector(c, m, m);} else {panic!("Matrix sizes do not match");}
}

        这种算法很慢,不仅因为三个嵌套循环,还因为按列通过而不是按行的内部循环遍历对于 CPU 缓存命中率来说是可怕的Bb.at(k, j)

三、换位以获得更好的性能

        转置朴素方法允许 B 上的乘法迭代在行而不是列上运行,将矩阵 B 的乘法步幅重新组织为更有利于缓存的格式。从而变成A x BA x B^t

         它涉及一个新的矩阵分配(无论如何,在这个实现中)和一个完整的矩阵迭代(一个 O(n²) 操作,更准确地说,这种方法是 O(n³) + O(n²))——我将进一步展示它的性能有多好。它如下所示:

fn multiply_transpose (A: Matrix, B: Matrix):C = new Matrix(A.num_rows, B.num_cols)// Construct transpose; requires allocation and iteration through BB’ = B.transpose()for i in 0 to A.num_rows:for j in 0 to B'.num_rows:sum = 0;for k in 0 to A.num_cols:// Sequential access of B'[j, k] is much faster than B[k, j]sum += A[i, k] * B'[j, k]C[i, j] = sumreturn C 

四、次立方:斯特拉森算法的工作原理

        要了解 Strassen 算法的工作原理(此处为 Rust 代码),首先考虑矩阵如何用象限表示。要概念化它的外观:

        在朴素算法中使用此象限模型,结果矩阵 C 的四个象限中的每一个都是两个子矩阵乘积的总和,总共产生 8 次乘法。

        考虑到这八个乘法,每个乘法都在一个块矩阵上运行,其行和列跨度约为 A 和 B 大小的一半,复杂性相同:

        斯特拉森算法定义了由这些象限组成的七个中间块矩阵

        仅通过 7 次乘法而不是 8 次乘法计算。这些乘法可以是递归斯特拉森乘法,并可用于组成最终矩阵:

由此产生的亚立方复杂度:

五、排比

        中间矩阵 M1 的计算 ...M7 是一个令人尴尬的并行问题,因此也很容易检测算法的并发变体(一旦你开始理解 Rust 关于闭包的规则)。

/*** Execute a recursive strassen multiplication of the given vectors, * from a thread contained within the provided thread pool.*/
fn 
_par_run_strassen (a: Vec<f64>, b: Vec<f64>, m: usize, pool: &ThreadPool) -> Arc<Mutex<Option<Matrix>>> {let m1: Arc<Mutex<Option<Matrix>>> = Arc::new(Mutex::new(None));let m1_clone = Arc::clone(&m1);pool.execute(move|| { // Recurse with non-parallel algorithm once we're // in a working threadlet result = mult_strassen(&mut Matrix::with_vector(a, m, m),&mut Matrix::with_vector(b, m, m));*m1_clone.lock().unwrap() = Some(result);});return m1;
}

六、标杆

        我编写了一些快速的基准测试代码,该代码在不断增加的矩阵维度范围内运行四种算法中的每一种进行几次试验,并报告每种算法的平均时间。

~/code/strassen ~>> ./strassen --lower 75 --upper 100 --factor 50 --trials 2running 50 groups of 2 trials with bounds between [75->3750, 100->5000]x    y    nxn      naive       transpose  strassen   par_strassen
75   100  7500     0.00ms      0.00ms     1.00ms     0.00ms
150  200  30000    6.50ms      4.00ms     4.00ms     1.00ms
225  300  67500    12.50ms     9.00ms     8.50ms     2.50ms
300  400  120000   26.50ms     22.00ms    18.00ms    5.50ms
[...]
3600 4800 17280000 131445.00ms 53683.50ms 21210.50ms 5660.00ms
3675 4900 18007500 141419.00ms 58530.00ms 28291.50ms 6811.00ms
3750 5000 18750000 154941.00ms 60990.00ms 26132.00ms 6613.00ms

        然后,我通过以下方式可视化结果:pyplot

        此图显示了矩阵从 7.5k 元素 () 到大约 19 万 () 的乘法时间。你可以看到朴素算法在计算上变得不切实际的速度有多快,在高端需要两分半钟。N x M = 75 x 100N x M = 3750 x 5000

        相比之下,Strassen 算法的扩展更平滑,并行算法计算两个 19M 个元素的矩阵的结果,而朴素算法只处理 3.6M 个元素所花费的时间。

        对我来说最有趣的是算法的性能。如前所述,缓存性能的改进(以牺牲完整矩阵副本为代价)在这些结果中得到了清楚地证明 - 即使使用与该方法渐近等效的算法也是如此。transposenaive

七、分析和性能优化

        这个文档是理解 Rust 性能基础知识的绝佳资源。在 Mac OS 上启动并运行仪器进行分析是微不足道的,这要归功于货运仪器的 Rust 指南。这是调查分配行为、CPU 热点和其他事情的绝佳工具。

在此过程中发生了一些变化:

  • Strassen 代码通过分而治之策略递归调用自己,但是一旦矩阵达到足够小的大小,其高常数因子使其比一般矩阵算法慢。我发现这个点是大约 64 的行宽或列宽;通过提高吞吐量提高几个因素来增加此阈值2
  • 斯特拉森算法要求矩阵填充到最接近的指数 2;减少这种情况以懒惰地确保矩阵只有偶数行和列 通过减少昂贵的大分配,将吞吐量提高了大约两倍
  • 将小矩阵回退算法从 更改为 导致大约 20% 的改进naivetranspose
  • 添加和添加到 Cargo.toml 发布构建标志大约提高了 5%。有趣的是,性能持续恶化codegen-units = 1lto = "thin"lto = “true”
  • 一丝不苟地删除所有可能的副本大约提高了~10%Vec
  • 提供一些提示并删除随机访问查找中的向量边界检查,又提高了大约 20%#[inline]
    /*** Returns the element at (i, j). Unsafe.*/#[inline]pub fn at (&self, i: usize, j: usize) -> f64 {unsafe {return *self.elements.get_unchecked(i * self.cols + j);}}

参考资料:

迈克·克维特
线性代数
算法

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/101299.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

【LLM数据篇】预训练数据集+指令生成sft数据集

note 在《Aligning Large Language Models with Human: A Survey》综述中对LLM数据分类为典型的人工标注数据、self-instruct数据集等优秀的开源sft数据集&#xff1a;alpaca_data、belle、千言数据集、firefly、moss-003-sft-data多轮对话数据集等 文章目录 note构造指令实例…

【Linux】网络层协议:IP

我们必须接受批评&#xff0c;因为它可以帮助我们走出自恋的幻象&#xff0c;不至于长久在道德和智识上自我陶醉&#xff0c;在自恋中走向毁灭&#xff0c;事实上我们远比自己想象的更伪善和幽暗。 文章目录 一、IP和TCP之间的关系&#xff08;提供策略 和 提供能力&#xff09…

QT 基本对话框

包括&#xff1a; 1.标准文件对话框 dialog.h #ifndef DIALOG_H #define DIALOG_H#include <QDialog> #include <QTextCodec> #include <QLabel> #include <QLineEdit> #include <QPushButton> #include <QGridLayout> #include <QFr…

通过DBeaver 给Postgre SQL表 设置主键自增

1.创建表 CREATE TABLE public.company ( id int4 NOT NULL , name text NOT NULL, age int4 NOT NULL, address bpchar(50) NULL, salary float4 NULL, join_date date NULL, CONSTRAINT company_pkey PRIMARY KEY (id) ); 2.插入数据&#xff08;不传入id&#xff…

【Leetcode】104.二叉树的最大深度

一、题目 1、题目描述 给定一个二叉树 root ,返回其最大深度。 二叉树的 最大深度 是指从根节点到最远叶子节点的最长路径上的节点数。 示例1: 输入:root = [3,9,20,null,null,15,7] 输出:3示例2: 输入:root = [1,null,2] 输出:2提示: 树中节点的数量在 [0, 104…

并查集及其简单应用

文章目录 一.并查集二.并查集的实现三.并查集的基本应用 一.并查集 并查集的逻辑结构:由多颗不相连通的多叉树构成的森林(一个这样的多叉树就是森林的一个连通分量) 并查集的元素(树节点)用0~9的整数表示,并查集可以表示如下: 并查集的物理存储结构:并查集一般采用顺序结构实…

【Redis】——Redis基础的数据结构以及应用场景

什么是redis数据库 Redis 是一种基于内存的数据库&#xff0c;对数据的读写操作都是在内存中完成&#xff0c;因此读写速度非常快&#xff0c;常用于缓存&#xff0c;消息队列、分布式锁等场景。&#xff0c;Redis 还支持 事务 、持久化、Lua 脚本、多种集群方案&#xff08;主…

Centos7卸载|安装JDK1.8|Xshell7批量控制多个终端

一: 使用yum安装的好处是较为方便|环境变量自动配置完成。 1.1: 执行下面的命令,检查是否已安装了jdk # 查看当前是否安装了JDK&#xff0c; [rootwww ~]# rpm -qa |grep java [rootwww ~]# rpm -qa |grep jdk [rootwww ~]# rpm -qa |grep gcj [rootwww ~]# rpm -qa | grep -…

数据结构——二叉搜索树(附带C++实现版本)

文章目录 二叉搜索树概念 二叉树的实际应用二叉树模拟实现存储结构二叉搜索树构成二叉搜索树的查找插入操作中序遍历二叉树的删除循环(利用左子树最右节点&#xff09;递归(利用右子树根节点) 二叉树拷贝二叉树资源的销毁 二叉树实现完整代码总结 二叉搜索树 概念 二叉搜索树…

PHPStudy 安装tp8 php8.2.9

一、PhpStudy升级PHP版本&#xff0c;安装PHP8.2操作步骤 1.1、官网下载最新的php版本 打开Windows版的官网下载&#xff0c;地址&#xff1a;https://windows.php.net/download/ 页面上有不同的PHP版本&#xff0c;这里我们下载的是64位nts版的PHP8.2.9。 1.2、解压下载的文…

2023.8 - java - 泛型

泛型问题的引出&#xff1a; jdk 1.5 引出泛型 // package 泛型; public class index {public static void main (String[] args){test t new test();t.setContent("aaa");int a (int) t.getContent();System.out.println(a);} }class test{Object content;publi…

国内常见的几款可视化Web组态软件

组态软件是一种用于控制和监控各种设备的软件&#xff0c;也是指在自动控制系统监控层一级的软件平台和开发环境。这类软件实际上也是一种通过灵活的组态方式&#xff0c;为用户提供快速构建工业自动控制系统监控功能的、通用层次的软件工具。通常用于工业控制&#xff0c;自动…

[SWPUCTF 2022 新生赛]ez_ez_php

这段代码是一个简单的PHP文件处理脚本。让我们逐行进行分析&#xff1a; error_reporting(0); - 这行代码设置了错误报告的级别为0&#xff0c;意味着不显示任何错误。 if (isset($_GET[file])) { - 这行代码检查是否存在一个名为"file"的GET参数。 if ( substr($_…

无涯教程-Perl - wantarray函数

描述 如果当前正在执行的函数的context正在寻找列表值,则此函数返回true。在标量context中返回false。 语法 以下是此函数的简单语法- wantarray返回值 如果没有context,则此函数返回undef&#xff1b;如果lvalue需要标量,则该函数返回0。 例 以下是显示其基本用法的示例…

记录Taro巨坑,找不到sass类型定义文件

问题 taronutuisassts项目里tsconfig.json一直报红报错。 找不到“sass”的类型定义文件。 程序包含该文件是因为: 隐式类型库 “sass” 的入口点 其实正常人想的肯定是装上 types/sass试试。开始我试过了&#xff0c;没用。删了依赖重装也没用。后面在issue中找到答案了 解决…

Mybatis的SqlSource SqlNode BoundSql

学习链接 MyBatis SqlSource解析 【Mybatis】Mybatis源码之SqlSource#getBoundSql获取预编译SQL Mybatis中SqlSource解析流程详解 Mybatis TypeHandler解析 图解 Mybatis的SqlSource&SqlNode - processon DynamicSqlSource public class DynamicSqlSource implement…

05-微信小程序常用组件-表单组件

05-微信小程序常用组件-表单组件 文章目录 表单组件button 按钮案例代码 form 表单案例代码 image 图片支持长按识别的码案例代码 微信小程序包含了六大组件&#xff1a; 视图容器、 基础内容、 导航、 表单、 互动和 导航。这些组件可以通过WXML和WXSS进行布局和样式设…

800V高压电驱动系统架构分析

需要电驱竞品样件请联&#xff1a;shbinzer &#xff08;拆车邦&#xff09; 过去一年是新能源汽车市场爆发的一年&#xff0c;据中汽协数据&#xff0c;2021年新能源汽车销售352万辆&#xff0c;同比大幅增长157.5%。新能源汽车技术发展迅速&#xff0c;畅销车辆在动力性能…

【Python原创设计】基于Python Flask的全国气象数据采集及可视化系统-附下载方式以及项目参考论文,原创项目其他均为抄袭

基于Python Flask的全国气象数据采集及可视化系统 一、项目简介二、项目技术三、项目功能四、运行截图五、分类说明六、实现代码七、数据库结构八、源码下载 一、项目简介 本项目是一个基于Web技术的实时气象数据可视化系统。通过爬取中国天气网的各个城市气象数据&#xff0c…

Nginx高可用集群

目录 一.简介二.案例1.实现思路2.配置文件修改3.实现效果故障转移机制 一.简介 以提高应用系统的可靠性&#xff0c;尽可能地减少中断时间为目标&#xff0c;确保服务的连续性&#xff0c;达到高可用的容错效果。例如“故障切换”、“双机热备”、“多机热备”等都属于高可用集…