算是给自己整理一下吧。主要是MF和FM的对比。

Toy Matrix#

假设有那么个用户数为6，物品数为5的大小为$6\times 5$的评分矩阵：

$$R=\begin{bmatrix}
4 & 0 & 5 & 1 & 0 \
3 & 5 & 1 & 0 & 0 \
2 & 5 & 0 & 0 & 2 \
0 & 0 & 0 & 5 & 4 \
1 & 5 & 5 & 0 & 0 \
0 & 0 & 4 & 4 & 3
\end{bmatrix}$$

空缺的评分用0填充。

评分预测#

MF模型#

这个是相对比较简单的，将MF的两个矩阵算一下点积，就能得到所有用户对所有物品的评分$\hat R$。

FM模型#

FM对输入数据$\bm x$（假设是行向量）的表示一般可以分为三部分：代表用户ID的one-hot向量$\bm{x}_u\in\mathbb{R}^6$，代表物品ID的one-hot向量$\bm{x}_i\in\mathbb{R}^5$，以及其他与该评分有关的上下文（context）特征向量$\bm{x}_c\in\mathbb{R}^c$（$c$为特征维度数）。即：$\bm x = [\bm x_u, \bm x_i, \bm x_c]$。

如果context是跟item挂钩的，即每个item下的所有评分都共享同一个context特征的话，事情就好办多了（至少这次跑的实验是这样的）。因此，FM的预测公式可以分解为与item相关以及与user相关两部分，跟MF模型相似。

原本计算FM的实现如下，$X$为行向量，$\bm w$为列向量，而$V$为$k\times(11+c)$的交叉项权重矩阵。

1

y_hat = w0 + X * w - sum((X .^ 2) * (V' .^ 2), 2) / 2 + sum((X * V') .^ 2, 2) / 2;

由于用户ID和物品ID都是one-hot向量，而context是跟item挂钩的，可以简化成：

1

y_hat = w0 + w(u) + w(i) + Xc(i) * wc - sum((Vu(u)' .^ 2) + (Vi(i)' .^ 2) + (Xc(i) .^ 2) * (Vc' .^ 2), 2) / 2 + sum((Vu(u)' + Vi(i)' + Xc(i) * Vc') .^ 2, 2) / 2

其中$u$和$i$是每条评分的用户和物品的ID。

然后跟MF一样，拆成与item挂钩的一部分：

1
2
3

linear_i = w_item + ctx * w_ctx; % [n_item, 1]
cross_i_pos = V_item' + ctx * V_ctx'; % [n_item, k]
cross_i_neg = V_item' .^ 2 + (ctx .^ 2) * (V_ctx' .^ 2); % [n_item, k]

ctx为n_item$\times c$的context特征矩阵。

以及跟user挂钩的另一部分：

1
2
3

linear_u = w_user; % [n_user, 1]
cross_u_pos = V_user'; % [n_user, k]
cross_u_neg = V_user' .^ 2; % [n_user, k]

最后对每个用户$u$的，预测公式如下：

1

y_hat = w0 + linear_u(u) + linear_i + sum(cross_u_pos(u) + cross_i_pos, 2)/2 - sum(cross_u_neg(u) + cross_i_neg, 2)/2; % [n_item, 1]

对用户没评过分的item，把itemID对应的位置置1，然后补上这个item的context特征进行预测就行了。

指标计算#

Precision@N#

N指的是推荐列表的长度。对推荐系统而言，前N个item的预测标签为1，而在后面的均为0。而对于ground truth标签，用户为其评过分则为1，否则为0。

For个example，假如$N=3$，推荐系统为第一个user生成的推荐列表为1，3，5，4，2。

因此TP=2，FP=1，Precision@3=TP/(TP+FP)=2/3=0.667。

实际上，计算可以简化为Pre@N=TP/N，TP为true positive的item数。代表的含义是前N个item列表中出现用户评过分的item的概率。

Recall@N#

跟上面类似，FN指的是出现在第N个item后面，用户有评过分的item数，如item 4。

因此，FN=1，Recall@3=TP/(TP+FN)=2/3=0.667。

实际上，计算也可以简化为Rec@N=TP/Np，Np为当前用户评过分的item个数。代表的含义是用户评过分的item出现在前N个item列表中的概率。这个数值是随N递增而单调递增的。

Average Precision#

这个指标综合评估Precision和Recall的质量，跟N无关。顾名思义，AP是对Precision求均值，而具体做法则是以recall为横坐标，precision为纵坐标，对precision求均值。人懒，原理就不解释了。

Normalized Discounted Cumulative Gain@N (NDCG@N)#

好像是这么拼来着吧？平时都直接喊的嗯滴希鸡。

推荐系统按照所有item预测的评分分数做个倒序排序，得到一个item列表：

Cumulative Gain，既然是cumulative，那自然就少不了sum嘛，CG@N=sum(2^r-1)，r按照生成item列表的顺序来。有些实现版本则是CG@N=sum(r)，没有取指数。

同样假设$N=3$，CG@3=(2^4-1)+(2^5-1)=46。

Discounted Cumulative Gain，多了个折损因子，DCG@N=CG@N/d。DCG@3=(2^4-1)/log2(2) + (2^5-1)/log2(3)=15+31/1.585=34.5588。

而Normalized则需要计算理想的DCG，即Ideal DCG（IDCG），它是按照评分倒序排的：

根据上表算DCG得到IDCG@3=DCG@3=(2^5-1)/log2(2) + (2^4-1)/log2(3) + (2^1-1)/log2(4)=31+15/1.585+1/2=40.9639。

最后的NDCG@N=DCG@N/IDCG@N=34.5588/40.9639=0.8436。

而某些论文中（这里就不点名是谁啦），NDCG的计算只考虑了用户评过分的物品，而未评过分的物品则不参与计算，因此实际计算NDCG时的列表如下（即去掉评分为0的item）：

计算IDCG时，是否去掉评分为0的item其实不影响结果，毕竟它们全排到最后面了嘛，求sum的时候都是0。

最后算出的DCG@3=15+31/1.585+1/2=35.0588。即NDCG@3=35.0588/40.9639=0.8558。用这种计算方法得出来的NDCG是偏高的，不过嘛，作者开心就好。

仓库地址： https://github.com/qhgz2013/shared_matrix

Bug该修的修好了，内存泄漏，多次free啥的操作基本上是没有了，跑小的数据集已经有挺好的效果的了：

RTX ON 共享内存 ON：
独立内存使用情况：

共享内存使用情况（在/dev/shm中）：

共享内存 OFF：

如果要跑大点的数据集，单用matlab自带的parallel.pool.Constant，就算是128G的内存也out of memory了，现在把15G的数据直接扔共享内存里，还算跑得动：

代码改动也不大，以前是：

1
2
3
4

x_const = parallel.pool.Constant(x);
parfor i = 1:n
    x_val = x_const.Value(:, i);
end

现在则是：

1
2
3
4
5
6
7
8

x_host = shared_matrix_host(x);
parfor i = 1:n
    x_dev = x_host.attach();
    x_val = x_dev.get_data();
    x_val = x_val(:, i);
    x_dev.detach();
end
x_host.detach();

好像是繁琐了一点，但是后面优化一下加个直接套个struct一次性创建多个变量的话应该简单很多，比如：

1
2
3
4
5
6
7
8
9

shared_memory_host = create_shmat_host(x, a, b, c);
parfor i = 1:n
    shared_memory_dev = attach_shmat_host(shared_memory_host);
    x_val = shared_memory_dev.x(:, i);
    a = shared_memory_dev.a;
    b = shared_memory_dev.b;
    detach_shmat_dev(shared_memory_dev);
end
detach_shmat_host(shared_memory_host);

这个就留给以后再说了。

参考：https://www.mathworks.com/matlabcentral/answers/454296-can-i-use-microsoft-visual-studio-2019-with-matlab-r2019a-or-r2018b

附件：MATLAB VS2019 Support.zip

步骤：

1. 将附件里的两个xml放到...\MATLAB\R2018b\bin\win64\mexopts下
2. 在注册表\HKEY_LOCAL_MACHINE\SOFTWARE\WOW6432Node\Microsoft\VisualStudio\SxS\VS7（找不到就新建）里添加字符串值：
   名称：16.0，
   数据：C:\Program Files (x86)\Microsoft Visual Studio\2019\Professional\（VS2019的目录）。

完毕，重启matlab即可。

目录#

共享内存的操作#

Windows下可以用Windows的API进行操作，而Linux下则可以通过POSIX的API进行操作，虽然过程有点区别，但终究是大同小异。

创建&数据写入#

为了简单测试，这里就直接用Windows版Python自带的mmap当实验：

1
2
3
4
5
6
7
8
9
10
11
12
13

>>> import numpy as np
>>> a = np.random.randn(32*4096,4096)  # 先占用4G内存
>>> a[0,0]
-0.5271477716515403
>>> a[-1,-1]
1.1194140920933267
>>> b = a.tobytes()  # 额外的4G内存，peak 8G
>>> del a  # 释放4G内存
>>> import mmap
>>> m = mmap.mmap(0, 32*4096*4096*8, 'Local\\MyShMem1')
>>> m.write(b)  # 额外的4G内存，peak 8G
>>> del b  # 释放4G内存
>>> 

这段代码会常占4G的内存，顶峰的内存会达到数据量的2倍。（在写入共享内存的时候）

注：Linux版Python不支持这种创建共享内存的方式。

数据读取#

数据读取是用matlab的mex写的：

shmem_win_data_read.c：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

#include <mex.h>
#include <stdio.h>
#include <windows.h>
#pragma comment(lib, "user32.lib")

TCHAR szName[] = TEXT("Local\\MyShMem1");
void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    //mxUint64 ptr = mxGetUint64s(prhs[0])[0];
    mxArray* pArr = mxCreateDoubleMatrix(32*4096, 4096, 0);
    HANDLE hMapFile = OpenFileMapping(FILE_MAP_ALL_ACCESS, FALSE, szName); // param: R/W, do not inherit the name, name of mapping object
    if (hMapFile == NULL) {
        printf("Failed to open file mapping: %d\n", GetLastError());
        return;
    }
    double* b = (double*)MapViewOfFile(hMapFile, FILE_MAP_ALL_ACCESS, 0, 0, 32*4096*4096*sizeof(double)); // Last: bufsize
    if (b == NULL) {
        printf("Failed to map view of file: %d\n", GetLastError());
        return;
    }
    for(int i = 0; i < 5; i++)
        printf("%lf\n", b[i]);
    double* originalArr = mxGetPr(pArr);
    mxFree(originalArr);
    printf("Original ptr: %lld\n", originalArr);
    mxSetPr(pArr, b);
    printf("Current ptr: %lld\n", mxGetPr(pArr));
    plhs[0] = pArr;
    if (nlhs > 1) {
        // handle
        mwSize d[] = { 1 };
        pArr = mxCreateNumericArray(1, d, mxUINT64_CLASS, mxREAL);
        *(unsigned long long*)mxGetData(pArr) = (unsigned long long)hMapFile;
        plhs[1] = pArr;
    }
}

然后编译，运行：

1
2
3
4
5

mex 'shmem_win_data_read.c' -g
a = cell(1);
[a{1}, h] = shmem_win_data_read;
b = a{1};
% 用变量b进行正常读写

返回的第一个参数是数据，至于为什么要用cell，就是为了规避matlab自带的GC机制。一旦将数组的指针利用mxSetPr设置为一个非通过调用mxMalloc得到的值（如c语言自带的malloc），GC机制尝试对它进行free的时候，后果自己可以想象。

返回的第二个参数是共享内存的handle，需要在使用完后释放掉。

cell的内存机制后面再说，注意的是，在调用mxSetPr之前，一定要用mxGetPr和mxFree把最初创建返回数组pArr时所申请到的内存释放掉，否则就会发生内存泄漏。

可以看出，数据是完全一致的。

资源释放#

一旦用完，就该妥善处理后事了（笑）。你大可以试试clear直接清掉，看看matlab会不会崩。

shmem_win_data_free.c：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

#include <mex.h>
#include <stdio.h>
#include <windows.h>
#pragma comment(lib, "user32.lib")
#include <matrix.h>
#include <stdlib.h>

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    mxArray* cell_arr = mxGetCell(prhs[0], 0);
    double* new_ptr = (double *)mxMalloc(sizeof(double) * 5);
    double* original_ptr = mxGetPr(cell_arr);
    new_ptr[0] = 1.2;
    new_ptr[1] = 0.0;
    new_ptr[2] = 0.0;
    new_ptr[3] = 0.0;
    printf("new ptr: %lld\n", new_ptr);
    printf("original ptr: %lld\n", original_ptr);
    mxSetPr(cell_arr, new_ptr);
    printf("current ptr: %lld\n", mxGetPr(cell_arr));
    mxSetM(cell_arr, 2);
    mxSetN(cell_arr, 2);
    UnmapViewOfFile(original_ptr);
    if (nrhs > 1) {
        HANDLE hMapFile = *(HANDLE*)mxGetData(prhs[1]);
        CloseHandle(hMapFile);
    }
}

为了简单实验起见，这里就直接用把内存改回matlab的托管内存了（用mxMalloc申请的内存），然后手动把之前改的内存释放掉。

1
2

mex 'shmem_win_data_free.c' -g
shmem_win_data_free(a, h);  % 这里就用a作为参数传入

合理利用matlab的“特性”，a{1}和b都会同时修改为一个新的数组，接着就可以直接clear掉，不需要再担心什么了。

完整使用#

1
2
3
4
5
6
7
8
9

v = zeros(1, 4096);
parfor x = 1:4096
    a = cell(1);
    [a{1}, h] = shmem_win_data_read;
    b = a{1};
    c = sum(b(:, x));
    v(x) = c;
    shmem_win_data_free(a, h);
end

这时候大可安心看系统的内存使用情况，妈妈再也不用担心内存占用啦。

番外：直接使用mxSetPr的效果#

mxSetPr在mex中一般用法效果如下：

1
2
3
4
5
6
7
8
9
10
11
12

#include <mex.h>
#include <stdio.h>

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    const mxArray* arg = prhs[0];
    double* ptr = mxGetPr(arg);
    printf("Original dataArray PTR: 0x%p\n", ptr);
    double* new_ptr = (double*)mxMalloc(sizeof(double) * mxGetM(arg) * mxGetN(arg));
    printf("Allocated dataArray PTR: 0x%p\n", new_ptr);
    mxSetPr(arg, new_ptr);
    printf("New dataArray PTR: 0x%p\n", mxGetPr(arg));
}

上面代码的作用：直接申请一块新的内存，将这块内存的地址直接赋值到matlab（作为输入参数传入mex函数）的矩阵数组中。运行结果如下：

可以看到，在mex函数体内对数组进行的任何修改是无法直接影响到输入参数的。

番外2：必须用cell中的元素作为输出变量#

还有一个关键点，就是必须用a{1}当输出变量，如：

1
2
3
4
5
6
7

a{1} = randn(5);
b = a{1};
% 通过b进行数据读取

do_cleanup_within_cell(a);
% 上面的函数能够同时作用于a（输入变量）和b（通过a引用赋值）
% 执行结果：a = 1×1 cell, a{1} = [], b = []

而把一般变量作为输出，再套用cell，则只能影响到cell里面的元素，不会影响到最先赋值的变量

1
2
3
4
5
6
7
8

t = randn(5); % 先赋值到一个新变量
a{1} = t; % 再通过引用赋值给cell
b = a{1};
% 通过t或者b进行数据读取

do_cleanup_within_cell(a);
% 上面的函数仅能作用于a（输入变量）和b（通过a引用赋值），而最原始的t则不受影响
% 执行结果：t = 5×5 double, a = 1×1 cell, a{1} = [], b = []

其中do_cleanup_within_cell可为如下（直接把输入数组清空为[]）：

1
2
3
4
5
6
7
8

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    mxArray* cell_arr = mxGetCell(prhs[0], 0);
    double* original_ptr = mxGetPr(cell_arr);
    printf("original ptr: 0x%p\n", original_ptr);
    mxSetM(cell_arr, 0);
    mxSetN(cell_arr, 0);
    mxSetPr(cell_arr, NULL);
}

反正个人猜测这种现象肯定跟matlab的GC是脱不了钩的，由于没有官方的内存介绍文档，所以也不深究了。写代码讲究一件事，那就是能用就行。

数据类型#

除了double、float以及(u)int8/16/32/64外，还有complex、sparse、cell、struct和object等数据类型，一个完整的check如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

#include <mex.h>
#include <matrix.h>
#include <stdio.h>

#define CHECK_FUNC(x) printf(#x": %d\n", x(prhs[0]))

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
    // 1. mxArray attributes
    CHECK_FUNC(mxIsNumeric);
    CHECK_FUNC(mxIsComplex);
    // 2. create, query, and access data types
    // 2.2 noncomplex float
    CHECK_FUNC(mxIsScalar);
    CHECK_FUNC(mxIsDouble);
    CHECK_FUNC(mxIsSingle);
    // 2.3 noncomplex integer
    CHECK_FUNC(mxIsInt8);
    CHECK_FUNC(mxIsUint8);
    CHECK_FUNC(mxIsInt16);
    CHECK_FUNC(mxIsUint16);
    CHECK_FUNC(mxIsInt32);
    CHECK_FUNC(mxIsUint32);
    CHECK_FUNC(mxIsInt64);
    CHECK_FUNC(mxIsUint64);
    // 2.6 sparse
    CHECK_FUNC(mxIsSparse);
    // 2.8 character
    CHECK_FUNC(mxIsChar);
    // 2.9 logical
    CHECK_FUNC(mxIsLogical);
    // 2.10 object
    // mxIsClass(mxArray* arr, char* className)
    // 2.11 structure
    CHECK_FUNC(mxIsStruct);
    // 2.12 cell
    CHECK_FUNC(mxIsCell);
}

接下来简单地验证完Linux的共享内存之后，就可以开始造（已经重复）的轮子了。

Linux版的碎碎念#

手头上只有Win版的R2018b和Linux版的R2017b服务器，想着版本差距不算大，直接改改代码，把WinAPI的Named shared memory改成POSIX的shm_open大概就行了。

把代码扔过去，编译，运行，哦豁完蛋，报了个在Windows下未曾出现过的崩溃错误：

突然告诉我还有个16Byte的header？喂不是吧大哥，我只把数据扔共享内存上了你告诉我header？

好吧，既然你想读，那就让你读呗。我把共享内存创建大一点，把数据扔后面一点，前面留16Byte的零，总算可以吧。然后：

你塔喵还要查header是吧？好好，我把mxGetPr往前16个字节的数据一并复制过去总行了吧。

嗯，这次好像真行了，而且能够正常退出没崩。

至于header是什么，说实话自己也没弄懂。随便弄了几个矩阵，往前挪16个字节print一下，结果如下：

131072*4096大小的double矩阵的header：
00 00 00 00 01 00 00 00 CE FA ED FE 20 00 10 00

7*65537大小的double矩阵：
40 00 38 00 00 00 00 00 CE FA ED FE 20 00 20 00

7*65538大小的double矩阵：
70 00 38 00 00 00 00 00 CE FA ED FE 20 00 10 00

7*65538大小的int32/uint32/single矩阵：
40 00 1C 00 00 00 00 00 CE FA ED FE 20 00 10 00

前面8字节明显跟数组的大小相关，那么后面8字节呢，跟type也不相关啊……特别是对sparse、complex这种数据类型也不变，而且倒数第二个0x10偶尔会变成0x20，也毫无规律可言。

目前只能推测跟matlab的内部GC有关（因为改一下清空变量的时候就会崩掉）。

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

目录#

引言#

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

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：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

#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的数组指针地址会发生变化。

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

1. zeros创建一个数组，matlab为其分配新内存，将引用数记为1，全部置0后返回这个数组
2. 创建变量a，将a中的数据指针的地址赋值为上面的数组地址，增加引用数（现在是2）
3. 清理调用zeros函数时创建的所有数组，包括减少zeros创建并返回的数组的一个引用数（现在是1），这条清理规则如果熟悉c++对象生命周期的话应该不难理解
4. randn创建一个数组，matlab为其分配新内存，将引用数记为1，赋值完成后返回这个数组
5. a变量的数据指针重新指向新的数组地址，新的数组的引用数增加为2，同时减少原本旧的数组的引用数（变为0）
6. 清理调用randn函数时创建的所有数组，这时将新数组的引用数减少为1
7. 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的时候，则会写入到原数组地址中。

栗子：

1
2
3
4

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.