上手实践ICCV2013的JDA（Joint Distribution Adaptation）方法

1、数据获取
　　选择surf特征文件作为算法的输入。surf特征文件可以从网盘上下载。下载到的文件主要包含4个.mat文件：amazon_SURF_L10.mat，Caltech10_SURF_L10.mat，dslr_SURF_L10.mat，webcam_SURF_L10.mat，它们对应 4 个不同的领域。彼此之间两两一组，就是一个迁移学习任务。每个数据文件包含两个部分：fts为800维的特征，labels为对应的标注。在测试中，选择由Caltech10_SURF_L10.mat作为源域，由amazon_SURF_L10.mat作为目标域。
　　Office+Caltech10数据集：Office是视觉迁移学习的主流基准数据集，包含3个对象领域Amazon(在线电商图片)、Webcam(网络摄像头拍摄的低解析度图片)、DSLR(单反相机拍摄的高解析度图片)，共有 4,652张图片 31 个类别标签。Caltech-256是对象识别的基准数据集，包括1个对象领域Caltech，共有 30,607张图片256个类别标签。对每张图片抽取SURF特征，并向量化为800维的直方图表征，所有直方图向量都进行减均值除方差的归一化处理，直方图码表由K均值聚类算法在Amazon子集上生成。具体共有4个领域C(Caltech-256), A(Amazon), W(Webcam) 和 D(DSLR)，从中随机选取2个不同的领域作为辅助领域和目标领域，则可构造4×3=12个跨领域视觉对象识别任务，如A→D，A→C，· · ·，C→W。
　　对surf文件夹中的4个数据文件使用如下代码进行简单归一化。具体方法是对数据进行加载，将最后的数据存入 Xs，Ys，Xt，Yt 这四个变量中。这四个变量分别对应源域的特征和标注、以及目标域的特征和标注。代码如下：

load("Caltech10_SURF_L10.mat");  % 源域
fts = fts ./ repmat(sum(fts,2),1,size(fts,2));  % fts是这个mat文件的变量，大小为N*M
Xs = zscore(fts,1); clear fts;  % z-score标准化
Ys = labels; clear labels;
save Caltech10_zscore_SURF_L10.mat Xs Ys;  % 保存到当前工作目录

load("amazon_SURF_L10.mat");  % 目标域
fts = fts ./ repmat(sum(fts,2),1,size(fts,2));  % fts是这个mat文件的变量，大小为N*M
Xt = zscore(fts,1); clear fts;  % z-score标准化
Yt = labels; clear labels;
save amazon_zscore_SURF_L10.mat Xt Yt;  % 保存到当前工作目录

这里只写出归一化待会要测试的Caltech10_SURF_L10.mat和amazon_SURF_L10.mat，在matlab工作区中建立Xs，Ys，Xt，Yt变量，若想测试其他数据集，可自行更改代码。

2、算法精炼
　　JDA主要进行边缘分布和条件分布的自适应。通过整理化简，JDA 最终的求解目标是：
$$(\mathrm{X}\sum _{\mathrm{c}=0}^{\mathrm{C}}{\mathrm{M}}_{\mathrm{c}}{\mathrm{X}}^{\mathrm{T}}+\lambda \mathrm{I})\mathrm{A}=\mathrm{X}\mathrm{H}{\mathrm{X}}^{\mathrm{T}}\mathrm{A}\mathrm{\Phi }$$ 上述表达式可以通过 Matlab 自带的eigs()函数直接求解。$\mathrm{A}$就是要求解的变换矩阵。下面明确各个变量所指代的含义:

• $\mathrm{X}$：由源域和目标域数据共同构成的数据矩阵
• $\mathrm{C}$：总的类别个数。在我们的数据集中，$\mathrm{C}=10$
• ${\mathrm{M}}_{\mathrm{c}}$： MMD矩阵。当$c=0$时为全MMD矩阵；当$c>1$时对应为每个类别的矩阵。
• $\mathrm{I}$：单位矩阵
• $\lambda$：平衡参数，直接给出
• $\mathrm{H}$：中心矩阵，直接计算得出
• $\mathrm{\Phi }$：拉格朗日因子，不用理会，求解用不到

3、编写代码
参考 JDA 开源的代码，直接给出精炼后的源码(来源于王晋东)

function [ acc,acc_iter,A ] = MyJDA( X_src,Y_src,X_tar,Y_tar,options )

% This is the implementation of Joint Distribution Adaptation.
% Reference: Mingsheng Long et al. Transfer feature learning with joint distribution adaptation.ICCV 2013.

% Inputs:
%%% X_src          :     source feature matrix, ns * n_feature
%%% Y_src          :     source label vector, ns * 1
%%% X_tar          :     target feature matrix, nt * n_feature
%%% Y_tar          :     target label vector, nt * 1

% Outputs:
%%% acc            :     final accuracy using knn, float
%%% acc_iter       :     list of all accuracies during iterations
%%% A              :     final adaptation matrix, (ns + nt) * (ns + nt)

% Set options
lambda = options.lambda;           %% lambda for the regularization
dim = options.dim;                 %% dim is the dimension after adaptation, dim <= n_feature
kernel_type = options.kernel_type; %% kernel_type is the kernel name, choose from 'primal' | 'linear' | 'rbf'
gamma = options.gamma;             %% gamma is the bandwidth of rbf kernel,can be missed for other 
T = options.T;                     %% iteration number, T >= 1. T <= 10 is suffice

acc_iter = [];
Y_tar_pseudo = [];
% Iteration
for i = 1 : T
    [Z,A] = JDA_core(X_src,Y_src,X_tar,Y_tar_pseudo,options);
    % normalization for better classification performance
    Z = Z*diag(sparse(1./sqrt(sum(Z.^2))));
    Zs = Z(:,1:size(X_src,1));
    Zt = Z(:,size(X_src,1)+1:end);

    knn_model = fitcknn(Zs',Y_src,'NumNeighbors',1);
    Y_tar_pseudo = knn_model.predict(Zt');
    acc = length(find(Y_tar_pseudo==Y_tar))/length(Y_tar); 
    fprintf('JDA+NN=%0.4f\n',acc);
    acc_iter = [acc_iter;acc];
end

end

function [ Z,A ] = JDA_core(X_src,Y_src,X_tar,Y_tar_pseudo,options)

% Set options
lambda = options.lambda;           %% lambda for the regularization
dim = options.dim;                 %% dim is the dimension after adaptation, dim <= m
kernel_type = options.kernel_type; %% kernel_type is the kernel name, primal|linear|rbf
gamma = options.gamma;             %% gamma is the bandwidth of rbf kernel

% Construct MMD matrix
X = [X_src',X_tar'];
X = X*diag(sparse(1./sqrt(sum(X.^2))));
[m,n] = size(X);
ns = size(X_src,1);
nt = size(X_tar,1);
e = [1/ns*ones(ns,1);-1/nt*ones(nt,1)];
C = length(unique(Y_src));

% M0
M = e * e' * C;  %multiply C for better normalization

% Mc
N = 0;
if ~isempty(Y_tar_pseudo) && length(Y_tar_pseudo) == nt
    for c = reshape(unique(Y_src),1,C)
        e = zeros(n,1);
        e(Y_src==c) = 1 / length(find(Y_src==c));
        e(ns+find(Y_tar_pseudo==c)) = -1 / length(find(Y_tar_pseudo==c));
        e(isinf(e)) = 0;
        N = N + e*e';
    end
end

M = M + N;
M = M / norm(M,'fro');

% Centering matrix H
H = eye(n) - 1/n * ones(n,n);

% Calculation
if strcmp(kernel_type,'primal')
    [A,~] = eigs(X*M*X'+lambda*eye(m),X*H*X',dim,'SM');
    Z = A'*X;
else
    K = kernel_jda(kernel_type,X,[],gamma);
    [A,~] = eigs(K*M*K'+lambda*eye(n),K*H*K',dim,'SM');
    Z = A'*K;
end

end

% With Fast Computation of the RBF kernel matrix
% To speed up the computation, we exploit a decomposition of the Euclidean distance (norm)
%
% Inputs:
%       ker:    'linear','rbf','sam'
%       X:      data matrix (features * samples)
%       gamma:  bandwidth of the RBF/SAM kernel
% Output:
%       K: kernel matrix
%
% Gustavo Camps-Valls
% 2006(c)
% Jordi (jordi@uv.es), 2007
% 2007-11: if/then -> switch, and fixed RBF kernel
% Modified by Mingsheng Long
% 2013(c)
% Mingsheng Long (longmingsheng@gmail.com), 2013

function K = kernel_jda(ker,X,X2,gamma)

switch ker
    case 'linear'
        if isempty(X2)
            K = X'*X;
        else
            K = X'*X2;
        end

    case 'rbf'
        n1sq = sum(X.^2,1);
        n1 = size(X,2);

        if isempty(X2)
            D = (ones(n1,1)*n1sq)' + ones(n1,1)*n1sq -2*(X'*X);
        else
            n2sq = sum(X2.^2,1);
            n2 = size(X2,2);
            D = (ones(n2,1)*n1sq)' + ones(n1,1)*n2sq -2*X'*X2;
        end
        K = exp(-gamma*D); 

    case 'sam'
        if isempty(X2)
            D = X'*X;
        else
            D = X'*X2;
        end
        K = exp(-gamma*acos(D).^2);

    otherwise
        error(['Unsupported kernel ' ker])
end

end

JDA方法的Matlab实现
将 JDA 方法包装成函数 MyJDA。
函数共接受 5 个输入参数：

• ${\mathrm{X}}_{\mathrm{s}\mathrm{r}\mathrm{c}}$：源域的特征，大小为$n_s\times m$
• ${\mathrm{Y}}_{\mathrm{s}\mathrm{r}\mathrm{c}}$：源域的标注，大小为$n_s\times 1$
• ${\mathrm{X}}_{\mathrm{t}\mathrm{a}\mathrm{r}}$：目标域的特征，大小为$n_t\times m$
• ${\mathrm{Y}}_{\mathrm{t}\mathrm{a}\mathrm{r}}$：目标域的标注，大小为$n_t\times 1$
• options：参数结构体，它包含：
• $\lambda c$：平衡参数，可以自由给出
• $T$：算法迭代次数
• $dim$：算法最终选择将数据降到多少维
• $kerneltype$：选择的核类型，可以选择RBF，线性、或无核
• $\gamma$：如果选择RBF核，那么它的宽度为$\gamma$

函数的输出包含3项：

• $acc$：算法的精度
• $ac{c}_{iter}$：算法每次迭代的精度，是一个一维数据
• $A$：最终的变换矩阵

4、测试算法
使用如下代码对JDA算法进行测试(可以单独建立.m文件，也可以直接在命令行下顺序输入该代码，确保所有文件均在同一目录下)：

options.T = 10;
options.gamma = 2;
options.kernel_type = 'linear';
options.lambda = 1.0;
options.dim = 20;
[Acc,Acc_iter,A] = MyJDA(Xs,Ys,Xt,Yt,options);
disp(Acc);

这里我将以上测试代码保存为test.m文件，在命令行窗口中敲下test或者软件上点运行。得到结果显示如下：

5、小结
通过以上过程，我们使用Matlab代码对JDA方法进行了实验，完成了一个迁移学习任务。其他的非深度迁移学习方法，均可以参考上面的过程。值得庆幸的是，许多论文的作者都公布了他们的文章代码，以方便我们进行接下来的研究。读者可以从Github或者相关作者的网站上获取其他许多方法的代码。

多说几句

运行在matlab R2017a，点击获取上面过程所有的代码。
快速上手步骤：

1. 打开matlab，定位到此文件夹。
2. 运行uniform.m，归一化源和目标域。
3. 运行test.m，可在工作区看到结果。