建公司网站,个人网站整站源码下载,百度排名点击器,征求网站建设矩阵平方和的计算
矩阵平方和的定义
矩阵平方和的定义是对矩阵中的每一个元素进行平方#xff0c;然后求和。
对于一个矩阵 A A A#xff0c;其平方和定义为#xff1a; sum ∑ i 1 m ∑ j 1 n A ( i , j ) 2 \text{sum} \sum_{i1}^{m}\sum_{j1}^{n} A(i,j)^2 sumi1∑…
矩阵平方和的计算
矩阵平方和的定义
矩阵平方和的定义是对矩阵中的每一个元素进行平方然后求和。
对于一个矩阵 A A A其平方和定义为 sum ∑ i 1 m ∑ j 1 n A ( i , j ) 2 \text{sum} \sum_{i1}^{m}\sum_{j1}^{n} A(i,j)^2 sumi1∑mj1∑nA(i,j)2
这个平方和计算在某些机器学习的算法中、或者特殊的优化问题中都会涉及到。
矩阵平方和的计算方法
通常而言前面的公式就给出了计算的方法无论怎么做都会涉及到对矩阵中的每一个元素进行平方然后求和。 因此算法的时间复杂度是 O ( m ∗ n ) O(m*n) O(m∗n)其中 m m m是矩阵的行数 n n n是矩阵的列数当 A A A是方阵时 m n mn mn。 最终算法的效率取决于矩阵的表示形式和对矩阵元素的访问方式。 当然还有可能会有一些特殊的优化方法比如矩阵的特殊性质可以通过一些特殊的方法来计算这里不做考虑。
最后就是并行指令集的使用比如SIMD指令集这里也不进行相关的讨论。
下面就主要是对在Matlab中实现矩阵平方和的几种方法进行比较。其实比较的是Matlab中访问矩阵元素的方式的性能差异。
Matlab的若干实现
我们非常随意的实现了一些方法。
其中最重要的就是所谓的基线方法。这个方法通常选择一个最简单/直观的方法便于抓住算法中的核心要素。 然后作为对比的基准也能够对不同算法的性能进行比较。
这种研究办法是一种常见的研究方法。例如在优化算法中通常选择格子搜索或者随机搜索作为基线方法。 提出一种新的方法通常希望能够对比基线方法有更好的性能并且这种性能提升是显著的。 此外也会结合新算法与基线算法的执行过程来解释算法优势的原理这就是一般算法研究文章中非常重要的理论分析。
function sumRet bench_loop_row_column(A)% 按照行、列的方式进行循环求和基线算法% 内存O(1)% 时间O(m*n) ~ O(n^2)[m, n] size(A);sumRet 0;% 对行进行循环for i 1:m% 对列进行循环for j 1:nsumRet sumRet A(i, j) ^ 2;end% 列循环结束end% 行循环结束
end上面这个最为直观的算法就是对矩阵的每一行进行遍历然后对每一行的每一个元素进行平方然后求和。 这个作为基线是一个恰当的选择。
相应的我们就会想到先循环列再循环行来作为对基线算法的改进。
function sumRet bench_loop_column_row(A)% 按照列、行的方式进行循环求和第一次改进的算法% 充分考虑到矩阵的存储方式列优先% 内存O(1)% 时间O(m*n) ~ O(n^2)[m, n] size(A);sumRet 0;for i 1:nfor j 1:msumRet sumRet A(j, i) ^ 2;endendend考虑到Matlab能够直接索引列我们就可以考虑对每列进行循环直接对列进行操作。或者换成对行进行操作。 下面几个算法就是行、列操作的形式这里又有一个小小的差别就是使用更多内存来直接进行向量加法最后调用sum求和。 或者还可以在循环内部调用sum函数这样内存的使用会从 O ( n ) \mathcal{O}(n) O(n)降低到
function s bench_loop_column_sum(A)% 直接循环列采用向量化的方式访问每个列% 用一个向量来存储每一列和% 最终调用sum函数求和% 内存O(n)% 时间考虑到向量加法的时间复杂度是O(n)% 所以这里的时间复杂度依然是O(m*n) ~ O(n^2)N size(A, 2);v zeros(size(A, 1), 1);for i 1:Nv v A(:, i) .^ 2;ends sum(v);
endfunction s bench_loop_sum_column(A)% 直接循环列采用向量化的方式访问每个列% 直接对列调用sum函数求列和最终累加求和% 内存O(1)% 时间O(m*n) ~ O(n^2)这里同样认为.^ 的时间复杂度是O(n)s 0;N size(A, 2);for i 1:Ns s sum(A(:, i) .^ 2);endendfunction s bench_loop_row_sum(A)% 直接循环行采用向量化的方式访问每行% 用一个向量来存储行和% 最终调用sum函数求和% 内存O(n)% 时间考虑到向量加法的时间复杂度是O(n)% 所以这里的时间复杂度依然是O(m*n) ~ O(n^2)N size(A, 1);v zeros(1, size(A, 2));for i 1:Nv v A(i, :) .^ 2;ends sum(v);
endfunction s bench_loop_sum_row(A)% 直接循环行采用向量化的方式访问每个行% 直接对行调用sum函数求行和最终累加求和% 内存O(1)% 时间O(m*n) ~ O(n^2)这里同样认为.^ 的时间复杂度是O(n)s 0;N size(A, 1);for i 1:Ns s sum(A(i, :) .^ 2);endend接下来我们利用矩阵的向量访问方式来构造一个算法。
function s bench_loop_vec(A)% 直接将矩阵展开成一个向量循环求和% 内存O(1)% 时间O(m*n) ~ O(n^2)s 0;N numel(A);for i 1:Ns s A(i) .^ 2;endend此外sum函数也提供了两个算法变体一个是调用两次默认对行求和
function s bench_sum_sum(A)% 直接两次调用sum函数对矩阵进行求和s sum(sum(A .^ 2));
endfunction s bench_sum_all(A)% 直接调用一次sum函数对矩阵进行求和s sum(A .^ 2, all);
end还可以把sum与矩阵列向量访问形式结合起来构成一种算法。
function s bench_sum_vec(A)% 直接将矩阵展开成一个向量% 进行元素.^计算调用sum函数求和s sum(A(:) .^ 2);
end最后还是利用矩阵的向量展开方式就是两个其实一样的算法因为dot函数内部就是调用矩阵乘法。
function s bench_vec_dot(A)% 直接将矩阵展开成一个向量进行点积计算s dot(A(:), A(:));
endfunction s bench_matrix_mul(A)% 直接将矩阵展开成一个向量进行矩阵乘法计算s A(:). * A(:);
end这些算法都非常平常但是通过上面的联系也能提升我们对Matlab的矩阵操作的理解。 接下来就是对这些算法的性能进行比较。
性能比较
性能比较方法
时间性能比较一般可以直接利用timeit函数来完成这个函数会对一个函数进行多次调用然后求中位数。 这个函数还能包括输出变量的构造时间。
利用这个函数我们编了一个工具对不同规模的矩阵进行比较。 这里代码的注释非常清楚无需赘述。
%% Benchmark functions helper
function [n, result] bench_f_n(n, f)% 按照给定的参数对函数进行测试% n: 矩阵大小% f: 函数句柄% 返回值n, result% n: 矩阵大小的向量 result: 运行时间的向量argumentsn (:, :) {mustBePositive}f function_handleend% 保证n是一个向量并且是正整数n round(n(:));result zeros(1, numel(n));% 对每一个n进行测试直接采用for循环不使用arrayfunfor i 1:numel(n)A rand(n(i), n(i));result(i) timeit(()f(A));end% result arrayfun((x) timeit(()f(rand(x, x))), n);% 这样也是可以的但是对于此处上面的写法更加直观% 这就是采用for循环的地方% 为了代码的可读性并且对性能的影响不大因为这里的循环次数不多
end性能比较代码
我们把前面的算法代码和上面的性能比较代码放在一个benchobjs的目录中就构成一个namespace命名空间通过import benchobjs.*来导入这些函数。 编写如下的测试脚本
Matlab code 设置测试范围运行测试获取数据可视化测试结果
%% import functions in benchobjs
import benchobjs.*;%% setup problem
n round(logspace(3, 4, 10));
functions {bench_loop_row_column, bench_loop_column_row, ...bench_loop_column_sum, bench_loop_sum_column, ...bench_loop_row_sum, bench_loop_sum_row, ...bench_loop_vec, bench_sum_sum, ...bench_sum_all, bench_sum_vec, ...bench_vec_dot, bench_matrix_mul};
markers {o, x, , *, s, d, ^, v, , , p, h};
% 这里是常见的小技巧将函数的句柄存储在一个cell数组中
% 这样可以方便的对这些函数进行遍历
% 同时把线型标签也存储在一个cell数组中这样可以方便的对这些线型进行遍历%% calculation
results cell(numel(functions), 2);for i 1:numel(functions)[n, result] bench_f_n(n, functions{i});results{i, 1} n;results{i, 2} result;fprintf(%s: %s: %s\n, func2str(functions{i}), ...mat2str(n), mat2str(result));
endcmd sprintf(save compareMatrixSquareSum%d results, maxNumCompThreads);
eval(cmd);% 这里试图考虑到多线程的影响目前在12核心的及其上进行了不同线程数的测试
% 发现对于这个问题多线程并没有展现出过大的差别
%% Visualize
% 一般也会把原始数据画出来稍微看一下
% 为了确保数据的正确性并且可以对数据进行初步的分析
figure;
clf;for i 1:numel(functions)[n, result] results{i, :};plot(n, result, LineWidth, 2, Marker, markers{i});yscale(log);hold on;fprintf(%s: %s: %s\n, func2str(functions{i}), ...mat2str(n), mat2str(result));
endlegend(cellfun((f)cellref(split(func2str(f), .), 2), ...functions, UniformOutput, false), ...Location, BestOutSide, interpreter, none);
xlabel(Matrix size);
ylabel(Time (s));
grid on% exportgraphics(gcf, ../matlab-img/compareMatrixSquareSum-time.png, ...
% Resolution, 600);%% Visualize 2
% 针对基准的加速比这是最常见的基准测试结果的展示方式
figure;
clf;for i 2:numel(functions)[n, result] results{i, :};plot(n, results{1, 2} ./ result, ...LineWidth, 2, Marker, markers{i});hold on;fprintf(%s: %s: %s\n, func2str(functions{i}), ...mat2str(n), mat2str(result));
endlegend(cellfun((f)cellref(split(func2str(f), .), 2), ...functions(2:numel(functions)), UniformOutput, false), ...Location, BestOutSide, interpreter, none);
xlabel(Matrix size);
ylabel(Time (s));
grid on% exportgraphics(gcf, ../matlab-img/compareMatrixSquareSum-acc.png, ...
% Resolution, 600);%% Visualize 2
% 去掉基准和两个最好的函数这样可以进行更加细节的比较和分析
% 这里必要性不太大主要是为了展示这种方式
figure;
clf;for i 2:numel(functions) - 2[n, result] results{i, :};plot(n, results{1, 2} ./ result, ...LineWidth, 2, Marker, markers{i});hold on;fprintf(%s: %s: %s\n, func2str(functions{i}), ...mat2str(n), mat2str(result));
endlegend(cellfun((f)cellref(split(func2str(f), .), 2), ...functions(2:numel(functions) - 2), UniformOutput, false), ...Location, BestOutSide, interpreter, none);
xlabel(Matrix size);
ylabel(Time (s));
grid on% exportgraphics(gcf, ../matlab-img/compareMatrixSquareSum-acc-2.png,...
% Resolution, 600);性能比较结果
最终可以得到如下的性能比较结果。 首先不同算法的性能有一定差异。基本上算法分为两组也可以视为三组。
第一组是基线算法对每一行进行遍历然后对每一个元素进行平方然后求和。第二组是向量展开访问并调用矩阵乘法mtimes直接调用BLAS二进制库的两个算法性能相当。其他就是各种循环的组合以及sum函数的组合。
这个分组按照加速比来看更加明显。 向量展开矩阵乘法的算法在1000~10000的规模下性能均有显著提升从15倍到40倍。
除去基线算法和向量化算法其他算法的关系较为复杂但是也能通过进行列访问、部分向量化来获得几倍的性能提升。 这充分显示了列优先矩阵访问对与计算效率的影响。 总结
进行算法开发一定要按照基线算法、算法优化的思路来考虑。对算法的效率进行比较最好选择不同的规模来分析问题。加速比是一个很好的指标能够直观的看出算法的性能提升。
