Вычислительные алгоритмы построения регрессии

Напишите программу, реализующую алгоритм Грама-Шмидта (алгоритм 11.3 в разделе 11.3). Опробуйте свою программу на данных Longley и на полиномиальных приближениях, описанных в примере 11.3 в разделе 11.8.4. Сравните свои результаты с результатами в примерах 11.2 и 11.3.

Проверьте формулы для алгоритма обновления в разделе 11.7.

Проведите несколько экспериментов, чтобы сравнить время выполнения регрессии, выполненной методом Холецкого, с регрессией, выполненной с помощью SVD. Дают ли приведённые в тексте подсчёты числа операций разумное представление об относительной скорости этих двух методов?

Вычислите алгебраически решение системы линейных уравнений, соответствующей матрице

используя гауссово исключение с частичным выбором ведущего элемента и без него. Покажите, что если частичный выбор ведущего элемента не используется, решение $x_{3}$ вычисляется как

а при использовании частичного выбора ведущего элемента — как $1-2 \epsilon$. Прокомментируйте точность этих двух методов, если $\epsilon$ очень мало.

Покажите, что применение одного шага гауссова исключения к матрице $\mathbf{X}^{T} \mathbf{X}$ даёт матрицу

где $\widetilde{\mathbf{X}}$ — центрированная версия $\mathbf{X}$.

Проверьте (11.15) и, следовательно, покажите, что $\sqrt{\mathrm{RSS}}$ задаётся $(p+1, p+1)$-м элементом матрицы $\mathbf{R}_{A}$.

Если $\mathbf{R}^{T} \mathbf{R}$ — разложение Холецкого матрицы $\mathbf{X}^{T} \mathbf{X}$, покажите, что

Пусть $\mathbf{q}_{1}, \ldots , \mathbf{q}_{p}$ — векторы, построенные в алгоритме 11.3. Покажите, что $\mathbf{q}_{1}, \ldots , \mathbf{q}_{j}$ — ортонормированный базис для $\mathcal{C}\left(\mathbf{a}^{(1)}, \ldots , \mathbf{a}^{(j)}\right)$ при $j = 1,2, \ldots , p$.

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

Выведите формулы (11.23), (11.24) и (11.25).

Докажите, что в начале этапа $j$ алгоритма MGSA столбцы $\mathbf{a}^{(1)}, \ldots , \mathbf{a}^{(j)}$ ортогональны, ортогональны $\mathbf{a}^{(j+1)}, \ldots , \mathbf{a}^{(p)}$ и порождают первые $j$ столбцов исходной матрицы.

Покажите, что произведение ортогональных матриц ортогонально.

Покажите, что если $\mathbf{Q}_{p}, \mathbf{Q}_{n-p}$ и $\mathbf{r}_{1}$ определены как в (11.36), то предсказанные значения $\widehat{\mathbf{Y}}$ задаются формулой

Пусть $\mathbf{G}_{i h}$ задано формулой $(11.40)$ и $\widetilde{\mathbf{D}} = \operatorname {diag}\left(\tilde{d}_{1}, \ldots , \tilde{d}_{n}\right)$. Покажите, что существует диагональная матрица $\mathbf{D} = \operatorname {diag}\left(d_{1}, \ldots , d_{n}\right)$ и быстрый ротатор вида (11.45), такие что $\mathbf{G}_{i h} = \mathbf{D M} \tilde{\mathbf{D}}^{-1}$. Покажите, что

где $t = s / c$ и $r = \tilde{d}_{h} / \tilde{d}_{i}$.

Докажите, что число обусловленности $\mathbf{X}_{w}$ задаётся формулой (11.61).

Повторите вычисления из примера 11.3, используя функцию $f(t) = \cos 4 t$ вместо $\exp (\sin 4 t)$. Сохраняется ли при этом ранжирование методов? Пример кода на R приведён после упражнения 3 ниже.

Измените приведённый ниже код так, чтобы подсчитывать число операций с плавающей точкой (flops), требуемых для выполнения каждого алгоритма. Насколько точны формулы, приведённые в таблице 11.2? ####################################################### # MGSA function: calculates Q and R^{-1} of QR decomp mgsa<-function(A){ n<-dim(A)[1] p<-dim(A)[2] GG<-diag(p) for(i in 1:(p-1)){ a<-A[,i] denom<-sum(aa) G<-numeric(p) for(j in (i+1):p){ num<-sum(A[,j]a) G[j]<- -num/denom A[,j]<-A[,j]+aG[j] } GG<-GG + outer(GG[,i],G) } list(W=A,G=GG) } ####################################################### #Householder function: calculates R of QR decomp of A house<-function(A){ n<-dim(A)[1] p<-dim(A)[2] for(j in 1:p){ indices<-j:n # construct Householder vector u<-numeric(n) u[indices]<-A[indices,j] unorm1<-sqrt(sum(u[indices]^2)) unorm2<-sum(u[indices[-1]]^2) uu<-u[j] u[j]<-if(uu<0)uu-unorm1 else -unorm2/(uu+unorm1) gamma<-0.5(unorm2 + u[j]^2) # multiply by householder matrix A[indices,j]<-0 A[j,j]<-unorm1 if(j==p)return(A) for(l in (j+1):p){ k<-sum(A[indices,l]u[indices])/gamma A[indices,l]<-A[indices,l] - ku[indices] } } A } ####################################################### # Givens function: calculates R of QR decomp of A givens<-function(A){ n<-dim(A)[1] p<-dim(A)[2] for(j in 1:(p-1)){ indices<-j:p for(i in (j+1):n){ sqr<-sqrt(A[j,j]^2+A[i,j]^2) cc<-A[j,j]/sqr ss<-A[i,j]/sqr temp<-ccA[j,indices] + ssA[i,indices] A[i,indices]<- -ssA[j,indices] + ccA[i,indices] A[j,indices]<-temp } } A }

Предположим, что мы подгоняем регрессию методом наименьших квадратов к подмножеству $J$ из $n$ имеющихся наблюдений. Пусть $\mathbf{X}_{J}$ — подматрица матрицы $\mathbf{X}$, соответствующая наблюдениям из $J$, пусть $e_{i}$ — $i$-й остаток из этой подгонки, $i = 1, \ldots , n$, и пусть $h_{i j} = \mathbf{x}_{i}^{T}\left(\mathbf{X}_{J}^{T} \mathbf{X}_{J}\right)^{-1} \mathbf{x}_{j}$. Покажите, что если мы удаляем наблюдение $i \in J$ и добавляем наблюдение $j \notin J$, то изменение остаточной суммы квадратов равно

Используя результат упражнения 1 выше, напишите функцию на R, реализующую метод Хокинса для вычисления приближения к оценке LTS.

Напишите функцию на R, реализующую алгоритм Рупперта для вычисления приближения к оценке LTS. Вам потребуется обратиться к работе Ruppert [1992] за деталями алгоритма. Разработайте небольшую симуляцию для сравнения алгоритма Рупперта с алгоритмом Хокинса.