2021-02-04

Matlab中基于共享内存的矩阵

（挖坑，目前只是试验可行性，不一定100%能做出什么东西出来）

目录#

引言
matlab的内存机制
matlab并行计算的内存机制

引言#

无论跑什么代码也好，只要数据的访问是相互独立，没有任何依赖关系的，一般情况下都可以引入并行计算增加CPU/GPU的利用率。matlab也不例外，最直接傻瓜的：将for改成parfor，就可以利用matlab自带的并行计算工具提高代码的运行效率。

举个栗子，就用手头上推荐系统的NDCG计算，每个用户的NDCG是相互独立的，所以为每个用户计算NDCG这一个过程就可以通过改成parfor i = 1:n实现并行化，最后再将结果取平均。

function ndcg = metric_ndcg(R, R_hat, k)
n = size(R, 1);
ndcg = zeros(1, n);
invalid_users = zeros(1, n);
RT = R';
RT_hat = R_hat';
parfor i = 1:n
    Ri = RT(:, i); % [1, items]
    Ri_hat = RT_hat(:, i)';
    non_empty_entries = find(Ri ~= 0);
    Ri = Ri(non_empty_entries);
    Ri_hat = Ri_hat(non_empty_entries);
    [~, idx] = sort(Ri_hat, 'descend');
    dcg_Ri = Ri(idx);
    dcg = dcg_at_k(dcg_Ri, k, true);
    
    idcg_Ri = sort(Ri, 'descend');
    idcg = dcg_at_k(idcg_Ri, k, true);
    if idcg ~= 0
        ndcg(i) = dcg / idcg;
    else
        invalid_users(i) = 1;
    end % if
end % for i
invalid_users = sum(invalid_users);
if n ~= invalid_users
    ndcg = sum(ndcg, 2)' / (n - invalid_users);
else
    ndcg = nan;
end
end % function

function dcg = dcg_at_k(r, k, p)
n = min([k length(r)]);
d = log2(2:n+1)';
rn = r(1:n);
if p
    rn = power(2, rn) - 1;
end
dcg = sum(rn ./ d);
end % function

当然，跑小的数据集不成问题，只不过评分矩阵R一旦大起来，内存的使用就会直线飙升。为什么？多亏了matlab的沙雕并行计算的内存机制，下面一一细数它的罪恶。

matlab的内存机制#

自己接触matlab不算多，平时也就写写代码，能跑就行的那种。但是，偶尔为了优化一下代码，摸点底层机制还是必要的。

为了探究这个机制，需要用到matlab中的mex，它是matlab提供的能调用其他编程语言的一个工具。这里就用c语言作为栗子，毕竟有它就能够直接实现很多骚操作：

官方文档：https://www.mathworks.com/help/matlab/matlab_prog/memory-allocation.html

mex c：从入门到治疗高血压#

首先需要#include <mex.h>，matlab调用c语言的函数void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])也是不可少的，其重要性相当于平时的main()。nlhs和nrhs两个值代表输入和输出的参数个数，输入的参数都保存在prhs数组中，而输出的参数则需要自己创建，赋值到plhs数组中。

这个栗子只是简单输出一下第一个输入参数的内存位置，函数体内第一行的mxArray PTR代表目前作为参数的这个矩阵对象的内存地址，而第二行dataArray PTR则是用mxGetPr获取这个对象里面的实际数据的内存位置，第三行就简单输出第一个数据的值。

编译和调用过程如下：

这两个指针的含义可以用一段等效的c++代码直观反映出来，在这个栗子中，mxArray PTR相当于&a，而dataArray PTR相当于a.array：

#include <string.h>

class mxArray {
public:
    double* array; // 实际保存数据的指针
    mxArray (int size_1, int size_2) { // 构造函数
        // ...
    }
    ~mxArray () {
        // ...
    }
};

// 调用的mxGetPr的等效代码
double* mxGetPr(const mxArray* arr) {
    return arr->array;
}

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    // ...
}

int main() {
    mxArray a(32*4096, 4096); // a = zeros(32*4096, 4096)
    const mxArray* b[] = { &a }; mexFunction(0, NULL, 1, b); // double_ptr(a)
    return 0;
}

矩阵的存储方式#

与python中的numpy不同，matlab是按维度索引从左到右进行存储的，比如a=randn(3,3)，在内存上则以a(1,1) a(2,1) a(3,1) a(1,2) a(2,2) a(3,2) a(1,3) a(2,3) a(3,3)的顺序分布，所以提高CPU缓存命中率的读取方式是a(:,1)这类固定后面若干维度的范围读取形式。

矩阵的GC机制#

matlab的GC是以数组的引用数为依据的，分配内存地址的时候以及进行读写操作时，matlab会自动跟踪并调整引用数，当引用数为0（意味着已经没有变量用到这段数据）的时候，GC会自动释放掉这段内存。

常见的赋值方式#

首先声明一下，赋值方式那名称都是我瞎起的……

为了更能掌握其中的数组指针变化，我写了一个新的mex函数用来创建数组：

第一种是最常见的直接赋值，即将右边的数组直接赋值到变量上：

1 2	a = new_arr; % <-- a = new_arr; % <--

并且第二条赋值语句完成后，变量a的数组指针地址会发生变化。

这种赋值操作的具体顺序是（当然不是每步都有验证，也没确定对不对）：

zeros创建一个数组，matlab为其分配新内存，将引用数记为1，全部置0后返回这个数组
创建变量a，将a中的数据指针的地址赋值为上面的数组地址，增加引用数（现在是2）
清理调用zeros函数时创建的所有数组，包括减少zeros创建并返回的数组的一个引用数（现在是1），这条清理规则如果熟悉c++对象生命周期的话应该不难理解
randn创建一个数组，matlab为其分配新内存，将引用数记为1，赋值完成后返回这个数组
a变量的数据指针重新指向新的数组地址，新的数组的引用数增加为2，同时减少原本旧的数组的引用数（变为0）
清理调用randn函数时创建的所有数组，这时将新数组的引用数减少为1
GC检测到有引用数为0的数组，释放原来那个数组的内存

第二种也是很常见的引用赋值，即将右边的数组没做任何修改直接赋值到新的变量上：

1 2	a = randn(32*4096, 4096); b = a; % <--

赋值完成后，变量b的数组指针与变量a一致。（同时这段数组的引用数会自增）至于说数据修改之后怎么影响，这就是matlab的另外一个机制：写入复制（copy on write）的事了。

第三种，一般很少用，我叫它覆盖赋值，前提是等号两边的数组个数要一致。赋值后，右边的数组数据会复制到左边数组中，不会改变左边的数组指针的地址：

1 2	a = new_arr; a(:) = new_arr; % <--

写入复制机制#

当数组的引用数大于1时，只要对其中一个引用了该数组的变量进行写入，都会触发写入复制机制，另外创建一个数组保存这个修改，原数组引用数-1。而引用数只有1的时候，则会写入到原数组地址中。

栗子：

a = new_arr;
b = a; % 引用赋值
a(1) = 233; % <--
b(1) = 666; % <--

然而有一个很有趣的“特性”：mex不受这个机制的约束，也就是说，只要拿到地址，就可以绕过这个机制进行写入，如：

当然，还有更绝的：

matlab并行计算的内存机制#

举个栗子，假如用E5的CPU跑parfor，数据矩阵A占4G，matlab会创建12个并行的计算进程占满12个CPU核心。整个架构管它叫master-slave也好，server-client也罢，反正master/server会向slave/client发送进行计算所需要的所有数据，slave/client独立进行计算，并将结果返回到master/server中。

那么沙雕matlab会怎么操作呢？

把数据矩阵A再复制12遍。没错，是再复制12遍。一个4G，那么加起来13个就会占掉52G的内存。如果看到这里已经感到不适，该去医院检查高血压了。

首先来几个FAQ：

Q1：既然matlab能调用c/c++代码，那我在c/c++代码里面手动创建线程，再调用matlab，不用它parfor不就行了吗？
A1：No，当你在c/c++创建新线程再调用matlab的时候，它已经注定要crash了。

Q2：那我在每次循环中只处理A中的某一行/某一列不就行了吗？
A2：No，就算每次只处理一行/列也好，两行/列也罢，matlab终究会把整个数据矩阵A完完整整地复制一遍的。

Q3：So easy，那我自己造个轮子，将共享内存的地址赋值到matlab，再手动释放，避开matlab的主动复制不就行了吗？
A3：可以，但愿你能从每次操作内存引发的crash中吸取教训。

因此，相对于从matlab自带的进程间通讯机制下手而言，选择相对容易解决的内存机制下手是比较明智的，好好利用matlab的一些“特性”应该不难实现。

To be continued.

雪饼的狗窝

黑贞天下第一