lus*_*oog 5 c transpose slice multidimensional-array matrix-multiplication
我正在构建一套函数来处理多维数组数据结构,我希望能够定义数组的任意切片,这样我就可以实现两个任意矩阵(又称Tensors或nd数组)的广义内积.
一个APL论文我看了(老实说,我无法找到它-我已经看了这么多)定义的矩阵产品上左矩阵X,其尺寸A;B;C;D;E;F和右矩阵Y,尺寸G;H;I;J;K,其中F==G作为
Z <- X +.× Y
Z[A;B;C;D;E;H;I;J;K] <- +/ X[A;B;C;D;E;*] × Y[*;H;I;J;K]
其中+/是和的和,并且×将元素逐个元素应用于两个相同长度的向量.
所以我需要左边的"行"切片和右边的"列"切片.我当然可以使用转置,然后使用"行"切片来模拟"列"切片,但我宁愿更优雅地做.
维基百科关于切片的文章引出了关于涂料载体的存根,这似乎是我正在寻找的奇迹治疗,但是没有太多可以继续下去.
如何使用涂料矢量来实现任意切片？
(很久以后我注意到了一个有一些细节的数组的Stride.)
通过dope向量或描述符引用每个数组,可以实现通用数组切片(无论是否构建到语言中) - 包含第一个数组元素地址的记录,然后是每个索引的范围和相应的系数.索引公式.该技术还允许立即数组转置,索引反转,子采样等.对于像C这样的语言,其中索引总是从零开始,具有d索引的数组的掺杂矢量具有至少1 + 2d参数.
http://en.wikipedia.org/wiki/Array_slicing#Details
这是一个密集的段落,但它实际上都在那里.所以我们需要一个像这样的数据结构:
struct {
TYPE *data; //address of first array element
int rank; //number of dimensions
int *dims; //size of each dimension
int *weight; //corresponding coefficient in the indexing formula
};
哪里TYPE是什么元素类型,该领域的矩阵.为简单和具体,我们将使用int.出于我自己的目的,我已经设计了各种类型的编码到整数句柄中,所以我int的工作是YMMV.
typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *arr;
所有指针成员都可以在与结构本身相同的已分配内存块中分配位置(我们将称之为标题).但是通过替换早期使用的偏移和struct-hackery,算法的实现可以独立于块内(或没有)内的实际内存布局.
自包含数组对象的基本内存布局是
rank dims weight data
dims[0] dims[1] ... dims[rank-1]
weight[0] weight[1] ... weight[rank-1]
data[0] data[1] ... data[ product(dims)-1 ]
间接数组共享数据(整个数组或1个或多个行切片)将具有以下内存布局
rank dims weight data
dims[0] dims[1] ... dims[rank-1]
weight[0] weight[1] ... weight[rank-1]
//no data! it's somewhere else!
并且包含沿另一个轴的正交切片的间接数组将具有与基本间接数组相同的布局,但是具有适当修改的调光和重量.
索引(i0 i1 ... iN)的元素的访问公式是
a->data[ i0*a->weight[0] + i1*a->weight[1] + ...
+ iN*a->weight[N] ]
假设每个索引i [j]在[0 ... dims [j])之间.
在weight通常布局的行主阵列的向量中,每个元素是所有较低维度的乘积.
for (int i=0; i<rank; i++)
weight[i] = product(dims[i+1 .. rank-1]);
因此,对于3×4×5阵列,元数据将是
{ .rank=3, .dims=(int[]){3,4,5}, .weight=(int[]){4*5, 5, 1} }
或者对于7×6×5×4×3×2阵列,元数据将是
{ .rank=6, .dims={7,6,5,4,3,2}, .weight={720, 120, 24, 6, 2, 1} }
因此,要创建其中一个,我们需要前一个问题中的相同辅助函数来计算维度列表中的大小.
/* multiply together rank integers in dims array */
int productdims(int rank, int *dims){
int i,z=1;
for(i=0; i<rank; i++)
z *= dims[i];
return z;
}
分配,只需malloc足够的内存并设置指向适当位置的指针.
/* create array given rank and int[] dims */
arr arraya(int rank, int dims[]){
int datasz;
int i;
int x;
arr z;
datasz=productdims(rank,dims);
z=malloc(sizeof(struct arr)
+ (rank+rank+datasz)*sizeof(int));
z->rank = rank;
z->dims = z + 1;
z->weight = z->dims + rank;
z->data = z->weight + rank;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
return z;
}
使用前一个答案中的相同技巧,我们可以创建一个变量参数接口,以便于使用.
/* load rank integers from va_list into int[] dims */
void loaddimsv(int rank, int dims[], va_list ap){
int i;
for (i=0; i<rank; i++){
dims[i]=va_arg(ap,int);
}
}
/* create a new array with specified rank and dimensions */
arr (array)(int rank, ...){
va_list ap;
//int *dims=calloc(rank,sizeof(int));
int dims[rank];
int i;
int x;
arr z;
va_start(ap,rank);
loaddimsv(rank,dims,ap);
va_end(ap);
z = arraya(rank,dims);
//free(dims);
return z;
}
甚至通过使用ppnarg的强大功能计算其他参数来自动生成秩论证.
#define array(...) (array)(PP_NARG(__VA_ARGS__),__VA_ARGS__) /* create a new array with specified dimensions */
现在构建其中一个非常容易.
arr a = array(2,3,4); // create a dynamic [2][3][4] array
使用函数调用检索元素elema,将每个索引与相应的权重相乘,对它们求和,并为data指针编制索引.我们返回一个指向元素的指针,以便调用者可以读取或修改它.
/* access element of a indexed by int[] */
int *elema(arr a, int *ind){
int idx = 0;
int i;
for (i=0; i<a->rank; i++){
idx += ind[i] * a->weight[i];
}
return a->data + idx;
}
相同的varargs技巧可以使界面更容易elem(a,i,j,k).
要做一个切片,首先我们需要一种方法来指定要提取的尺寸和要折叠的尺寸.如果我们只需要从维度中选择单个索引或所有元素,那么该slice函数可以将rank int作为参数,并将-1解释为选择整个维度或0 .. dims i -1选择单个索引.
/* take a computed slice of a following spec[] instructions
if spec[i] >= 0 and spec[i] < a->rank, then spec[i] selects
that index from dimension i.
if spec[i] == -1, then spec[i] selects the entire dimension i.
*/
arr slicea(arr a, int spec[]){
int i,j;
int rank;
for (i=0,rank=0; i<a->rank; i++)
rank+=spec[i]==-1;
int dims[rank];
int weight[rank];
for (i=0,j=0; i<rank; i++,j++){
while (spec[j]!=-1) j++;
if (j>=a->rank) break;
dims[i] = a->dims[j];
weight[i] = a->weight[j];
}
arr z = casta(a->data, rank, dims);
memcpy(z->weight,weight,rank*sizeof(int));
for (j=0; j<a->rank; j++){
if (spec[j]!=-1)
z->data += spec[j] * a->weight[j];
}
return z;
}
收集与参数数组中的-1s对应的所有维度和权重,并用于创建新的数组头.所有参数> = 0都乘以它们的相关权重并添加到data指针,使指针偏移到正确的元素.
就数组访问公式而言,我们将其视为多项式.
offset = constant + sum_i=0,n( weight[i] * index[i] )
因此,对于我们选择单个元素(+所有较低维度)的任何维度,我们将现在常量索引的因子分解并将其添加到公式中的常量项(在我们的C表示中是data指针本身) .
辅助函数casta使用shared创建新的数组头data.slicea当然,改变权重值,但通过计算权重本身,casta成为更普遍可用的功能.它甚至可以用于构造一个动态数组结构,该结构直接在常规的C风格多维数组上运行,从而进行转换.
/* create an array header to access existing data in multidimensional layout */
arr casta(int *data, int rank, int dims[]){
int i,x;
arr z=malloc(sizeof(struct arr)
+ (rank+rank)*sizeof(int));
z->rank = rank;
z->dims = z + 1;
z->weight = z->dims + rank;
z->data = data;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}
return z;
}
涂料矢量也可用于实现转座.维度(和索引)的顺序可以更改.
请记住,这不像其他人那样正常的"转置".我们根本不重新排列数据.这更恰当地称为'dope-vector伪转置'.我们只需更改访问公式中的常量,重新排列多项式的系数,而不是更改数据.从实际意义上讲,这只是对交换性和加法相关性的应用.
因此,为了具体,假设数据从假设地址500开始顺序排列.
500: 0
501: 1
502: 2
503: 3
504: 4
505: 5
506: 6
如果a是等级2,变暗{1,7),权重(7,1),则索引的总和乘以添加到初始指针(500)的相关权重,产生每个元素的适当地址
a[0][0] == *(500+0*7+0*1)
a[0][1] == *(500+0*7+1*1)
a[0][2] == *(500+0*7+2*1)
a[0][3] == *(500+0*7+3*1)
a[0][4] == *(500+0*7+4*1)
a[0][5] == *(500+0*7+5*1)
a[0][6] == *(500+0*7+6*1)
因此,dope-vector伪转置重新排列权重和维度以匹配索引的新排序,但总和保持不变.初始指针保持不变.数据不会移动.
b[0][0] == *(500+0*1+0*7)
b[1][0] == *(500+1*1+0*7)
b[2][0] == *(500+2*1+0*7)
b[3][0] == *(500+3*1+0*7)
b[4][0] == *(500+4*1+0*7)
b[5][0] == *(500+5*1+0*7)
b[6][0] == *(500+6*1+0*7)
因此,通过使用通用切片或转置+"行" - 更复杂(更容易),可以实现广义内积.
首先,我们需要两个辅助函数,用于对产生矢量结果的两个矢量应用二元运算,并使用二元运算减少矢量,从而产生标量结果.
与上一个问题一样,我们将传入运算符,因此可以将许多不同的运算符使用相同的函数.对于这里的样式,我将运算符作为单个字符传递,因此已经存在从C运算符到函数运算符的间接映射.这是一个x-macro表.
#define OPERATORS(_) \
/* f F id */ \
_('+',+,0) \
_('*',*,1) \
_('=',==,1) \
/**/
#define binop(X,F,Y) (binop)(X,*#F,Y)
arr (binop)(arr x, char f, arr y); /* perform binary operation F upon corresponding elements of vectors X and Y */
表中的额外元素reduce用于null向量参数的函数.在这种情况下,reduce应该返回运算符的标识元素,0表示+,1表示*.
#define reduce(F,X) (reduce)(*#F,X)
int (reduce)(char f, arr a); /* perform binary operation F upon adjacent elements of vector X, right to left,
reducing vector to a single value */
因此binop,对操作符进行循环和切换.
/* perform binary operation F upon corresponding elements of vectors X and Y */
#define BINOP(f,F,id) case f: *elem(z,i) = *elem(x,i) F *elem(y,i); break;
arr (binop)(arr x, char f, arr y){
arr z=copy(x);
int n=x->dims[0];
int i;
for (i=0; i<n; i++){
switch(f){
OPERATORS(BINOP)
}
}
return z;
}
#undef BINOP
如果有足够的元素,则reduce函数会执行向后循环,如果有的话,将初始值设置为最后一个元素,并将初始值预设给运算符的identity元素.
/* perform binary operation F upon adjacent elements of vector X, right to left,
reducing vector to a single value */
#define REDID(f,F,id) case f: x = id; break;
#define REDOP(f,F,id) case f: x = *elem(a,i) F x; break;
int (reduce)(char f, arr a){
int n = a->dims[0];
int x;
int i;
switch(f){
OPERATORS(REDID)
}
if (n) {
x=*elem(a,n-1);
for (i=n-2;i>=0;i--){
switch(f){
OPERATORS(REDOP)
}
}
}
return x;
}
#undef REDID
#undef REDOP
使用这些工具,内部产品可以以更高级别的方式实施.
/* perform a (2D) matrix multiplication upon rows of x and columns of y
using operations F and G.
Z = X F.G Y
Z[i,j] = F/ X[i,*] G Y'[j,*]
more generally,
perform an inner product on arguments of compatible dimension.
Z = X[A;B;C;D;E;F] +.* Y[G;H;I;J;K] |(F = G)
Z[A;B;C;D;E;H;I;J;K] = +/ X[A;B;C;D;E;*] * Y[*;H;I;J;K]
*/
arr (matmul)(arr x, char f, char g, arr y){
int i,j;
arr xdims = cast(x->dims,1,x->rank);
arr ydims = cast(y->dims,1,y->rank);
xdims->dims[0]--;
ydims->dims[0]--;
ydims->data++;
arr z=arraya(x->rank+y->rank-2,catv(xdims,ydims)->data);
int datasz = productdims(z->rank,z->dims);
int k=y->dims[0];
arr xs = NULL;
arr ys = NULL;
for (i=0; i<datasz; i++){
int idx[x->rank+y->rank];
vector_index(i,z->dims,z->rank,idx);
int *xdex=idx;
int *ydex=idx+x->rank-1;
memmove(ydex+1,ydex,y->rank);
xdex[x->rank-1] = -1;
free(xs);
free(ys);
xs = slicea(x,xdex);
ys = slicea(y,ydex);
z->data[i] = (reduce)(f,(binop)(xs,g,ys));
}
free(xs);
free(ys);
free(xdims);
free(ydims);
return z;
}
此函数还使用cast提供varargs接口的函数casta.
/* create an array header to access existing data in multidimensional layout */
arr cast(int *data, int rank, ...){
va_list ap;
int dims[rank];
va_start(ap,rank);
loaddimsv(rank,dims,ap);
va_end(ap);
return casta(data, rank, dims);
}
它还用于vector_index将1D索引转换为索引的nD向量.
/* compute vector index list for ravel index ind */
int *vector_index(int ind, int *dims, int n, int *vec){
int i,t=ind, *z=vec;
for (i=0; i<n; i++){
z[n-1-i] = t % dims[n-1-i];
t /= dims[n-1-i];
}
return z;
}
github文件.其他支持函数也在github文件中.
此Q/A对是一个系列的,其出现在落实相关问题的一部分,印加了用C写的别人APL语言的解释:我如何与动态分配的任意维数组工作？,以及如何将C数学运算符(+ - */%)传递给函数result = mathfunc(x,+,y);? .其中一些材料以前发布在comp.lang.c中.comp.lang.apl中的更多背景信息.
| 归档时间： |
|
| 查看次数： |
925 次 |
| 最近记录： |