Метод наименьших углов (пример)

Материал из MachineLearning.

Перейти к: навигация, поиск

Метод наименьших углов (англ. least angle regression, LARS) - алгоритм отбора признаков в задачах линейной регрессии. Алгоритм предложили Бредли Эфрон, Тревор Хасти и Роберт Тибширани [1].

LARS решает следующую задачу. Пусть назначена линейная регрессионная модель. При большом количестве свободных переменных возникает проблема неустойчивого оценивания весов модели. LARS предлагает метод выбора такого набора свободных переменных, который имел бы наиболее значимую статистическую связь с зависимой переменной. Также LARS предлагает метод оценки весов.

Алгоритм LARS похож на алгоритм шаговой регрессии. Различие в том, что алгоритм LARS, вместо последовательного добавления свободных переменных, на каждом шаге изменяет их веса. Веса увеличиваются так, чтобы доставить наибольшую корреляцию с вектором регрессионных остатков. Основным достоинством LARS является то, что он выполняется за число шагов, не превышающее число свободных переменных. Частным случаем алгоритма LARS является алгоритм LASSO.

Содержание

Постановка задачи

Задана выборка - матрица X, столбцы которой соответствуют независимым переменным, а строки - элементам выборки и вектор \mathbf{y}, содержащий элементы зависимой переменной. Назначена линейная модель \mathbf{y}=X\mathbf{\beta}+\mathbf{\varepsilon}.

Обозначим множество столбцов матрицы X как \{\mathbf{x}_1,\ldots,\mathbf{x}_n\}. Пусть выполнены условия нормировки

\|\mathbf{y}\|_1=0,     \|\mathbf{x}_j\|_1=0,     \|\mathbf{x}_j\|_2^2=1,    где    j=1,\ldots,n.

Предполагается, что векторы \mathbf{x}_1,\ldots,\mathbf{x}_n линейно независимы.

Критерием качества модели назначена среднеквадратичная ошибка

S(\mathbf{\beta}) = \|\mathbf{y}-\mathbf{\mu}\|^2,     где   \mathbf{\mu} = X\mathbf{\beta}.

Требуется найти для каждого шага LARS вектор весов \mathbf{\beta} и набор столбцов, доставляющий наибольшую корреляцию векторов \mathbf{\mu} и \mathbf{y}.

Описание алгоритма

LARS последовательными шагами строит оценку коэффициентов \hat{\mathbf{\beta}}. На k-м шаге только k элементов вектора \hat{\mathbf{\beta}} отличны от нуля. Алгоритм последовательно вычисляет приближение вектора значений зависимой переменной

\mathbf{\mu}=X\mathbf{\beta}.

Для приближений используется вектор корреляций столбцов матрицы X с вектором остатков \mathbf{y}-\mathbf{\mu}

\mathbf{c}(\mathbf{\mu})=X^T(\mathbf{y}-\mathbf{\mu}).
Пример работы алгоритма в случае двух переменных. Пусть вектор  является проекцией вектора  на линейное подпространство . Назначим начальное приближение . Вектор регрессионных остатков  коррелирует с вектором  больше, чем с вектором . Первый шаг заключается в оценке . Скаляр  выбирается таким образом, что вектор остатков  делит пополам угол между векторами  и . Далее получаем значение , где  - нормированный вектор, делящий этот угол пополам. Так как мы рассматриваем случай двух переменных, то .
Пример работы алгоритма в случае двух переменных. Пусть вектор \bar{\mathbf{y}}_2 является проекцией вектора \mathbf{y} на линейное подпространство \text{span}(\mathbf{x}_1,\mathbf{x}_2). Назначим начальное приближение \hat{\mathbf{\mu}}_0= </p> \mathbf{0}. Вектор регрессионных остатков \bar{\mathbf{y}}_2-\hat{\mathbf{\mu}}_0 коррелирует с вектором \mathbf{x}_1 больше, чем с вектором \mathbf{x}_2. Первый шаг заключается в оценке \hat{\mathbf{\mu}}_1=\hat{\mathbf{\mu}}_0+\gamma_1\mathbf{x}_1. Скаляр \gamma_1 выбирается таким образом, что вектор остатков \bar{\mathbf{y}}_2-\hat{\mathbf{\mu}}_0 делит пополам угол между векторами \mathbf{x}_1 и \mathbf{x}_2. Далее получаем значение \hat{\mathbf{\mu}}_2=\hat{\mathbf{\mu}}_1+\gamma_2\mathbf{u}_2, где \mathbf{u}_2 - нормированный вектор, делящий этот угол пополам. Так как мы рассматриваем случай двух переменных, то \hat{\mathbf{\mu}}_2=\bar{\mathbf{y}}_2.

На k-м шаге новое значение приближения вектора значений зависимой переменной \mathbf{y} вычисляется как

\hat{\mathbf{\mu}}_k=\hat{\mathbf{\mu}}_{k-1}+\gamma_k\mathbf{u}_k.

Здесь \mathbf{u}_k - нормированный вектор, задающий биссектрису между добавляемым вектором — столбцом матрицы X и вектором регрссионных остатков.

Он вычисляется следующим образом. Пусть \mathcal{A}  - подмножество индексов \{1,\dots,n\} столбцов матрицы X. Это подмножество задает подматрицу

X_\mathcal{A} = [s_{j_1}\mathbf{x}_{j_1},\ldots,s_{j_{|\mathcal{A}|}}\mathbf{x}_{j_{|\mathcal{A}|}}], j\in\mathcal{A},

где множитель s\in\{+1, -1\} и |\mathcal{A}| - мощность множества \mathcal{A}. Обозначим ковариационную матрицу

\mathcal{G} = X_{\mathcal{A}}^TX_{\mathcal{A}}   и   A_{\mathcal{A}} = (\mathbf{1}_{\mathcal{A}}^T\mathcal{G}^{-1}_{\mathcal{A}}\mathbf{1}_{\mathcal{A}})^{-\frac{1}{2}},

где \mathbf{1}_{\mathcal{A}}  - вектор, состоящий из |\mathcal{A}| единиц.

Вычислим единичный вектор

\mathbf{u}_{\mathcal{A}} = X_{\mathcal{A}}\mathbf{w}_{\mathcal{A}}, </p>   где  \mathbf{w}_{\mathcal{A}} = A_{\mathcal{A}}\mathcal{G}_{\mathcal{A}}^{-1}\mathbf{1}_{\mathcal{A}}.

Вектор \mathbf{u} образует со столбцами матрицы X_{\mathcal{A}} одинаковые углы, меньшие \frac{\pi}{2}. Справедливы равенства

X'_{\mathcal{A}}\mathbf{u}_{\mathcal{A}} = A_{\mathcal{A}}\mathbf{1}_{\mathcal{A}}  и  \|\mathbf{u}_{\mathcal{A}}\|^{2} = 1.

Выполнение алгоритма

Назначим начальную оценку вектора значений зависимой переменной \hat{\mathbf{\mu}}=\mathbf{0}. Вычислим текущую оценку \hat{\mathbf{\mu}}_{\mathcal{A}} и вектор корреляций

\hat{\mathbf{c}} = X^T(\mathbf{y}-\hat{\mathbf{\mu}}_{\mathcal{A}}).

Найдем текущий набор индексов \mathcal{A}, соответствующих признакам с наибольшими абсолютными значениями корреляций

\hat C = \max_{j=1,\ldots,n}\{|\hat c_j|\}   и  \mathcal{A} = \{j:|\hat c_j| = \hat C\}.

Пусть s_j = \text{sign}(\hat c_j) для j\in\mathcal{A}. Построим матрицы X_{\mathcal{A}}, A_{\mathcal{A}}. Вычислим вектор \mathbf{u}_{\mathcal{A}} и вектор скалярных произведений

\mathbf{a}=X^T\mathbf{u}_{\mathcal{A}}.

Пересчитаем значение вектора \hat{\mathbf{\mu}}_{\mathcal{A}}

\hat{\mathbf{\mu}}_{\mathcal{A_{+}}} = \hat{\mathbf{\mu}}_{\mathcal{A}} + \hat\gamma\mathbf{u}_{\mathcal{A}},

где

\hat\gamma = {\min_{j\in\mathcal{A}^c}}^{+}\lbrace\frac{\hat{C}-\hat c_j}{A_{\mathcal{A}}-a_j},\frac{\hat C + \hat c_j}{A_{\mathcal{A}}+a_j}\rbrace. (*)

Минимум берется по всем положительным значениям аргументов для каждого j.

Алгоритм повторяется n раз.

Замечание

Величина \gamma в формуле (*) интерпретируется следующим образом. Определим

\mathbf{\mu}(\gamma) = {\mathbf{\mu}}_{\mathcal{A}}+\gamma\mathbf{u}_{\mathcal{A}}

при условии \gamma>0. Корреляция по добавляемому j-му признаку равна

c_j(\gamma)=\mathbf{x}_j^T(\mathbf{y}-\mathbf{\mu}(\gamma)) = c_{j} - \gamma a_{j}.

Для j\in\mathcal{A} получаем

|c_{j}(\gamma)|=C-\gamma A_{\mathcal{A}}.

Это означает, что все рассматриваемые на данном шаге максимальные абсолютные корреляции уменьшаются на одну и ту же величину. Из предыдущих двух соотношений мы видим, что при j\in \mathcal{A}^c корреляция c_j(\gamma) принимает наибольшее значение при

\gamma = \frac{C-c_{j}}{A_{\mathcal{A}}-a_j}.

Аналогично, корреляция -c_j(\gamma) принимает наибольшее значение при

\gamma = \frac{C+c_{j}}{A_{\mathcal{A}}+a_j}.

Таким образом, \gamma в выражении (*) - минимальная положительная величина, при которой новый индекс j может быть добавлен в набор \mathcal{A}.

Для иллюстрации работы алгоритма задумаем модель f = w_1+w_2\sqrt(x)+w_3x (более подробно см.). 

f = inline('w(1)*x.^0 + w(2)*x.^(1/2) + w(3)*x.^1', 'w','x'); % synthetic model
b = [1 -2 2]'; % parameters
x = [0.1:0.1:1]'; % independent variable
y = f(b,x); % dependent variable

Пусть имеется выборка, в которой вектор-столбцы матрицы X

X = \left( \begin{array}{ccccc} x^0_1 & \sqrt{x}_1 & x\sqrt{x}_1 & x_1 & x\log{x}_1 \\ \vdots& \vdots & \vdots &\vdots &\vdots \\ x^0_m & \sqrt{x}_m & x\sqrt{x}_m & x_1 & x\log{x}_m \\ \end{array} \right)

и вектор \mathbf{y}=(y_1,\ldots,y_m). 

% We use the following features to find the model we thought of
g = {'x.^0', 'x.^(1/2)', 'x.^1', 'x.^(3/2)', 'x.*log(x)'};
% make the sample set; columns are independent variables
X = eval([ '[' g{1} ' ' g{2} ' ' g{3} ' ' g{4} ' ' g{5} ']']);

Принята линейная модель \mathbf{y}=X\mathbf{\beta}. Требуется последовательно отобрать вектор-столбцы матрицы X, которые имеют наибольшую корреляцию с вектором остатков. Алгоритм LARS в данном примере сделает n=5 шагов. 

[n,p] = size(X);
muA = zeros(n,1);  % estimation of the dependent variable
beta = zeros(p,1); % estimation of parameters
betaLtst = beta';  % keep parameters in a storage
 
for i = 1:p
    % correlation coefficients between each feature (column of X) and vector of residuals
    c = X'*(y-muA); % note that columns of X are centered and normalized
    [C, A] = max(abs(c)); % find maximal value of correlation and corresponding index j of column in X
    % A = find(C == abs(c));
    % Aplus = find(C==c); % never used
    Sj  = sign(c(A)); % get sign of j-th correlation coefficient
    XA = X(:,A).*(ones(n,1)*Sj'); %
    % XA = X(:,A)*Sj; %
    G = XA'*XA; % norm of XA
    oA = ones(1,length(A)); % vector of ones
    AA =(oA*G^(-1)*oA')^(-0.5); % inverse matrix in the normal equation
 
    wA = AA*G^(-1)*oA'; % parameters to compute the unit bisector
    uA = XA*wA; % compute the unit bisector
    a = X'*uA; % product vector to compute new gamma
    if i<p % for all columns of X but the last
        M = [(C-c)./(AA-a);(C+c)./(AA+a)];
        M(find(M<=0)) = +Inf;
        gamma = min(M);
    else
        gamma  = C/AA;
    end
 
    muA = muA + gamma*uA; % make new approximation of the dependent variable
    beta(A) = beta(A) + gamma*wA.*Sj; % make new parameters
    betaLtst = [betaLtst; beta']; % store the parameters at k-th step
end

Результатом работы алгоритма является последовательность добавляемых признаков (веса которых отличны от нуля). 

   0.5739         0         0         0         0
   0.5739    0.0311         0         0         0
   0.5739    0.0311         0         0         0
   0.5739    0.0439         0         0         0
   0.5739    0.0439    0.2135         0         0
На графике показано изменение вектора параметров по шагам LARS. По оси абсцисс отложена сумма параметров, по оси ординат - их значения. Согласно исходной модели, только первые три параметра отличны от нуля.
На графике показано изменение вектора параметров \mathbf{\beta} по шагам LARS. По оси абсцисс отложена сумма параметров, по оси ординат - их значения. Согласно исходной модели, только первые три параметра отличны от нуля.
Ось абсцисс - свободная переменная, ось ординат - зависимая. Точками показаны исходные данные, соответствующие задуманной модели. Красная линия - модель LARS с параметрами , полученными с помощью LARS. Красная линия - модель для параметров, выбранных LARS и настроеных с помощью метода наименьших квадратов.
Ось абсцисс - свободная переменная, ось ординат - зависимая. Точками показаны исходные данные, соответствующие задуманной модели. Красная линия - модель LARS с параметрами \mathbf{\beta}, полученными с помощью LARS. Красная линия - модель для параметров, выбранных LARS и настроеных с помощью метода наименьших квадратов.

Литература

  • Tibshirani R. Regression shrinkage and Selection via the Lasso. // Journal of the Royal Statistical Society.Series B(Metodological). 1996  Vol. 32, \No 1, p.267-288.
  • Efron B., Hastie T., Johnstone J., Tibshirani R. Least Angle Regression. // The Annals of Statistics, 2004. Vol. 32, \No 2, p. 407-499.

См. также

Источник — «http://machinelearning.ru/wiki/index.php?title=%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%BD%D0%B0%D0%B8%D0%BC%D0%B5%D0%BD%D1%8C%D1%88%D0%B8%D1%85_%D1%83%D0%B3%D0%BB%D0%BE%D0%B2_%28%D0%BF%D1%80%D0%B8%D0%BC%D0%B5%D1%80%29»
Личные инструменты