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

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

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=0,$ nbsp; где 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=\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_\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}},$ где $\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_i}{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% 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 с параметрами $\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.

См. также

Полный текст данного примера
Регрессионный анализ
Регрессионная модель
Линейная регрессия (пример)
Порождение линейных регрессионных моделей (постановка задачи)
Лассо

Теги: Метод наименьших углов , Регрессионный анализ , модель LARS , алгоритм LASSO

Автор: Strijov | 23.04.2009 13:00 |

генератор прогнозов \| краткосрочный прогноз потребления электроэнергии \| финансовый прогноз \| прогнозирование в экономике \| онлайн прогнозирование
	© strijov.ru, 2008-2009 \| контакты: strijov@gmail.com \| карта сайта

Генератор прогнозов создает модель с помощью обучающих данных и набора элементарных функций. Результат работы генератора - прогноз и структура модели данных.