Линейная регрессия. Какой алгоритм использовать для решения метода наименьших квадратов — обратный, LU или ?

Я работаю над алгоритмом для выполнения линейной регрессии для одной или нескольких независимых переменных.

то есть: (если у меня есть m значения реального мира и в случае двух независимых переменных a и b)

C + D*a1 + E* b1 = y1

C + D*a2 + E* b2 = y2

...

C + D*am + E* bm = ym

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

Я буду использовать матричные обозначения, поэтому

введите здесь описание изображения

где Beta — это вектор [C, D, E], где эти значения будут наиболее подходящими линиями.

Вопрос Как лучше всего решить эту формулу? Должен ли я вычислять инверсию введите здесь описание изображения

или я должен использовать факторизацию/декомпозицию матрицы LU. Какова производительность каждого из них при большом объеме данных (т.е. большое значение m может быть порядка 10 ^ 8...)

ИЗМЕНИТЬ

Если ответ заключался в использовании разложения Холецкого или разложения QR, есть ли какие-либо подсказки реализации/простые библиотеки для использования. Пишу на C/C++.


person Saher Ahwal    schedule 21.10.2013    source источник


Ответы (2)


Ваша матрица X^TX должна иметь разложение Холецкого. Я бы изучил эту декомпозицию перед LU. Это быстрее: http://en.wikipedia.org/wiki/Cholesky_decomposition

person Tom Swifty    schedule 21.10.2013
comment
как насчет вычисления обратного? я должен избегать этого? -- Спасибо - person Saher Ahwal; 21.10.2013
comment
Обычно предпочтительнее избегать инверсии матрицы. Вы можете столкнуться с проблемами численной стабильности, и память может быть неэффективной для хранения инверсии. Придерживайтесь методов декомпозиции (они существуют не просто так). Методы Холецкого можно найти во многих пакетах (gsl и т. д.). - person Tom Swifty; 21.10.2013
comment
Спасибо. Я ценю ваш ответ - person Saher Ahwal; 21.10.2013
comment
Какой язык вы используете? - person Fernando; 21.10.2013
comment
@Saher, из-за размера вашей матрицы, вы также можете прочитать этот пост: stackoverflow.com/questions/13148052/ - person Tom Swifty; 21.10.2013
comment
Я использую C, а также буду использовать CUDA C на графическом процессоре. - person Saher Ahwal; 22.10.2013
comment
протестировали общую линейную регрессию в Python и использовали Cholesky. Спасибо - person Saher Ahwal; 24.10.2013

На ум приходят два простых подхода к решению плотной переопределенной системы Ax=b:

  1. Сформируйте A^T A x = A b, затем разложите по Холески A^T A = L L^T, затем выполните два обратных решения. Обычно это дает вам ответ с точностью до sqrt (машинный эпсилон).

  2. Вычислите факторизацию QR A = Q*R, где столбцы Q ортогональны, а R является квадратным и верхнетреугольным, используя что-то вроде исключения Хаусхолдера. Затем решите Rx = Q^T b вместо x обратной подстановкой. Это обычно дает вам ответ с точностью примерно до машинного эпсилон --- точность в два раза выше, чем у метода Холецкого, но это занимает примерно вдвое больше времени.

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

person tmyklebu    schedule 21.10.2013