Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
где /? — невырожденная верхняя треугольная матрица. '
Так как Q — ортогональная матрица, то из условия(23.18) следует ! *
RTu = 0, (23.22)
sTu + tTv - со > 0. (23.23)
Поскольку R не вырождена, равенство (23.22) означает, что и = 0. Поэтому (23.23) сводится к
Гги = со>0. (23.24)
Из (23.21) вытекает, что л-я компонента х„ вектора х сама является псевдорешением редуцированной задачи
Псевдообратной матрицей для столби** !HMtW» Bfp»atl tTl(tTt). Поэтому решение задачи (23.25) можно ВМЯШСЯТЬЯI
со
77
>0.
(23.26)
Лемма 23.17 доказана.
Ж
Можно считать, что алгоритм NNLS состоит из основного цикла (цикла А) и внутреннего цикла (цикла В). Цикл В составляют шаги 6-11; единственный вход в цикл В — через шаг 6, а единственный выход— через шаг 7.
Цикл А состоит из шагов 2—5 и цикла Д. Он начинается шагом 2; выход из цикла - через шаг 3.
На шаге 2 множество ^определяет положительные компоненты текущего вектора х. Компоненты х, индексированные множеством 26, в этот момент равны нулю.
Выбранный на шаге 4 индекс t указывает компоненту, не представленную пока в множестве &, которая в соответствии с леммой 23.17 будет положительна, если ее ввести в решение. На шаге 6 как раз и происходит ввод этой компоненты в пробное решение z. Если все прочие компоненты z, индексированные множеством 9* , остаются положительными, то на шаге 7 выполняется присваивание х: =z и возврат к началу цикла А. В рассматриваемом случае множество & расширяется, а 26 уменьшается эа счет переноса индекса t.
Во многих задачах попросту происходит повторение этой последовательности событий с добавлением очередной положительной компоненты при каждом выполнении цикла А; в конечном счете произойдет выход через шаг 3.
Если, однако, некоторая компонента, индекс которой находится в множестве Ьъ, стала неположительной в векторе z шага 6, то шаг 7 вынуждает алгоритм остаться в цикле В и перейти от х к х + a(z -х), 0 < а < 1. Множитель а выбирается как можно большим; нужно лишь сохранить неотрицательность нового вектора .v. Цикл В повторяется, пока в конце концов не произойдет выход через шаг 7.
Конечность цикла В можно установить, заметив, что все операции внутри этого цикла определены корректно, что при каждом выполнении шага 11 по крайней мере один индекс, а именно индекс, называемый здесь q, удаляется из множества & и что zr всегда положительно (это следует из леммы 23.17). Если я - число индексов в множестве & при входе в цикл В, то выход через шаг 7 должен произойти не позднее, чем через и - 1 повторений цикла. На практике выход из цикла ? обычно происходит при первом же достижении шага 7; при этом шаги 8-11 не выполняются вовсе.
Конечность цикла А можно установить, заметив, что норма невязкн
р(*)= \f~Ex В
строго уменьшается при каждом очередном достижении шага 2. Тем самым вектор х и ассоциированное с ним множество &-{i: xt > 0 } на шаге 2 каждый раз иные. Так как Фестъ подмножество в {1,...,//) и таких подмножеств конечное число, то цикл А должен закончиться после конечного числа итераций. На ряде малых тестовых задач было замечено, что обычно цикл А требует около л/2 итераций.
Модификация QR-pазложения матрицы И. Задача наименьших квадратов, решаемая на шаге 6, отличается от задачи, решавшейся при предыдущем выполнении этого шага, либо тем, что на шаге 5 в задачу был включен дополнительный столбец, либо тем, что один или несколько столбцов Е были удалены на шаге 11. Если сохранено ДО-разложение преды-
126 i
душей задачи, то для вычисления нового QR-разложения можно применить приемы модификации. Три метода модификации будут описаны в гл. 24.
Учет эффектов машинной арифметики. Когда шаг 6 выполняется сразу вслед за шагом 5, то вычисленное на шаге 6 значение компоненты zt теоретически должно быть положительно. Если жег, вследствие погрешностей округлений неположительно, то может произойти деление на нуль на шаге 8 или, возможно, будет неправильно вычислено а = 0 на шаге 9.
Этих неприятностей можно избежать, проверяя значение z, после шага 6 в тех случаях, когда вход в этот шаг был из шага 5. Если окажется, что z г < 0, это можно интерпретировать как указание, что значение w,, вычисленное на шаге 2 и участвовавшее в проверках шагов 3 и 4, следует скорее считать нулем, чем положительным числом. Поэтому можно положить wt: = 0 и вернуться к шагу 2. Результатом будет или выход через шаг 3, или присвоение t нового значения на шаге 4.
На шаге 11 всякое xit вычисленное значение которого отрицательно (что может произойти только вследствие округлений), должно рассматриваться как нулевое, а его индекс нужно переместить иэ множества & в множество 2?-
Проверка знаков zt,i е & , на шагах 7 и 8, по-видимому, не является критической. Возможные погрешности в классификации не должны иметь серьезных последствий.
Решение задачи LDP (23.3) можно получить путем подходящей нормировки вектора невязки соответствующей задачи NNLS (23.2). На этот метод решения задачи LDP (и на его обоснование) внимание авторов обратил Аллан Клайн.
Пусть даны m X л-матрица G, целые числа тили т-вектор //. Если система неравенств Gx> л совместна, то логической переменной *р в алгоритме будет присвоено значение TRUE и будет вычислен вектор v минимальной длины, удовлетворяющий этой системе. Если же неравенства несовместны, то алгоритм полагает ip = FALSE, и вектору х не приписывается никакого значения. Рабочие массивы, необходимые алгоритму, не включены в его список параметров.
