CC BY
18
УДК 519.237.5:519.24
МЕТОДЫ УСТОЙЧИВОГО ПОСТРОЕНИЯ ЛИНЕЙНЫХ МОДЕЛЕЙ НА ОСНОВЕ СПУСКА ПО УЗЛОВЫМ ПРЯМЫМ1
А. Н. Тырсин, А. А. Азарян
RESISTANT LINEAR MODEL FITTING METHODS BASED ON THE DESCENT THROUGH THE NODAL STRAIGHT LINES
А. N. Tyrsin, A. A. Azarian
Аннотация. Актуальность и цели. При построении линейных моделей во многих случаях приходится сталкиваться со стохастической неоднородностью экспериментальных данных. Это проявляется в нарушении условий теоремы Гаусса-Маркова, в частности, наблюдения могут быть засорены грубыми ошибками. В этих условиях оценивание параметров моделей требуется выполнять с помощью устойчивых методов. К их числу относят метод наименьших модулей и обобщенный метод наименьших модулей. Однако известные алгоритмы их реализации являются достаточно эффективными лишь для малых размерностей моделей и ограниченного объема выборок. Цель данного исследования - разработка эффективных вычислительных алгоритмов реализации данных методов, не имеющих ограничений на порядок моделей и объем экспериментальных данных. Материалы и методы. Реализация задач была достигнута за счет организации спуска к искомому решению по узловым прямым. Для снижения вычислительных затрат использована особенность узловых прямых - все расположенные на каждой такой прямой узловые точки являются пересечением набора гиперплоскостей, из которых отличной является только одна гиперплоскость. Результаты. Дан обзор известных методов реализации метода наименьших модулей и обобщенного метода наименьших модулей. Описаны алгоритмы спуска по узловым прямым, позволяющие существенно снизить вычислительные затраты при реализации этих методов. Доказано достижение минимума за конечное число шагов при реализации метода наименьших модулей с помощью спуска по узловым прямым. Реализация обобщенного метода наименьших модулей с помощью спуска по узловым прямым позволяет найти глобальный минимум или близкое к нему решение. Выводы. Предложены эффективные алгоритмы реализации метода наименьших модулей и обобщенного метода наименьших модулей при оценивании параметров линейных моделей, основанные на спуске по узловым прямым. Вычислительная сложность этих алгоритмов делает возможной их практическую реализацию для анализа экспериментальных данных и построения многомерных линейных моделей.
Ключевые слова: метод наименьших модулей, обобщенный метод, линейная модель, регрессия, алгоритм, узловая точка, узловая прямая, гиперплоскость.
Abstract. Background. When fitting linear models, in many cases one has to deal with stochastic inhomogeneity of the experimental data. This is manifested in violation of the assumptions of the Gauss-Markov theorem, in particular, observations can contain outliers. Under the circumstances the estimation of the parameters of models is required to perform using resistant methods. Among them are the least absolute deviations method and the generalized method of the least absolute deviations. However, the known algorithms for
1 Работа выполнена при финансовой поддержке гранта РФФИ, проект № 16-06-00048а.
their implementation are sufficiently effective only for small dimensions of models and a limited sample size. The purpose of this study is the development of effective computational algorithms for the implementation of these methods, which have no limitations on the order of the models and the amount of experimental data. Materials and methods. The implementation of the tasks was achieved by the descent to the solution through the nodal straight lines. To reduce computational costs, the feature of nodal straight lines is used - all nodes located on each such straight line are intersections of a set of hyperplanes, of which only one hyperplane is different. Results. A review of known methods for realizing the least absolute deviations method and the generalized method of the least absolute deviations is given. The algorithms of descent through the nodal straight lines are described, which allow to significantly reduce the computational costs when implementing these methods. The achievement of a minimum in a finite number of steps when implementing the method of the least absolute deviations based on descent through the nodal straight lines is proved. The implementation of the generalized method of the least absolute deviations by descent through the nodal straight lines allows us to find a global minimum or a solution close to it. Conclusions. Effective algorithms for realization of the least absolute deviations method and the generalized method of the least absolute deviations when estimating the parameters of linear models based on descent through the nodal straight lines are proposed. The computational complexity of these algorithms makes possible their practical implementation for the analysis of experimental data and the construction of multidimensional linear models.
Key words, least absolute deviations method, generalized method, linear model, regression, algorithm, nodal point, nodal straight line, hyperplane.
Введение
При исследовании стохастических систем часто используют линейные модели. К этому классу моделей относят линейные регрессионные зависимости, динамические модели авторегрессии и распределенного лага, а также линейные стохастические модели временных рядов (модели авторегрессии, авторегрессии - скользящего среднего).
Модель линейной регрессии имеет вид [1]
где
X {xij )nxm
( 1 1
y = Xa + £,
x Л
(x л
v xn J
(1)
заданная матрица;
(a \
a =
V am J
e Rm
( * л *2
Уп
fc ^
e Rn
£ =
c
Vе n J
e Rn — векторы коэффициентов регрес-
сии, измерений и случайных ошибок, соответственно.
Модель авторегрессии и распределенного лага, обозначаемая как ЛОЬ( р, д1,..., дт), записывается в виде [2]
у,=ao+Z aky—k+Z Z j i—j+ci •
к=1 l=1 j=1
(2)
Если в (2) р = 0 , то получим модель распределенного лага ОЬ^,...,qm)
т q¡
У* = а0 + ЕЕЬ1Х1,,-j + е,- • (3)
1=1 j=l
При т = 0 , а0 = 0 имеем модель авторегрессии АЯ(р)
р
У, = Е акУ-к + е, , (4)
к =1
которая при введении линейной комбинации лагов случайной компоненты будет моделью авторегрессии-скользящего среднего АЯМА(р, q):
р q
у1=Е акУ-к+Е bj • (5)
к=1 1=1
Построение конкретной математической модели по имеющимся наблюдениям реализуется с помощью статистических методов оценки ее параметров. Многие задачи, связанные с обработкой статистических данных, решаются в предположении существования достаточной информации об изучаемых объектах, процессах, явлениях и о свойствах, действующих на них возмущений. Для широкого класса задач разработаны методы эффективного оценивания неизвестных параметров с использованием классических методов максимального правдоподобия [3, 4]. В частности, в предположении, что случайные ошибки нормально распределены, методом максимального правдоподобия является метод наименьших квадратов (МНК) [1].
На основе МНК создана целостная система статистической обработки. С учетом простоты реализации он является наиболее распространенным статистическим методом построения линейных моделей. Его использование для оценивания параметров модели (1) требует выполнения ряда предпосылок, называемых условиями Гаусса-Маркова [5]. При их выполнении МНК-оценки параметров модели (1) являются состоятельными и несмещенными. Кроме того, если случайные ошибки имеют нормальный закон распределения, то МНК-оценки модели (1) становятся эффективными. У динамических моделей (2)-(5) условие некоррелированности значений случайной компоненты для различных лагов не выполняется, поэтому МНК-оценки остаются состоятельными и несмещенными.
Использование МНК при нарушении условий Гаусса-Маркова может привести к значительным ошибкам при оценивании параметров, а в случае присутствия в измерениях больших выбросов - даже к несостоятельности оценок [5]. В этом случае устойчивость оценок может быть достигнута за счет использования иных методов, самым популярным из которых является метод наименьших модулей (МНМ) [6-9]. МНМ для задачи (1) имеет вид
Q(a) = ЕУ - (а,^ . (6)
. . 1 1 aeR
,=1
Если выбросы в целом имеют симметричный вид, то для оценивания моделей (1), (3), не использующих в качестве входных переменных лагов вы-
ходной величины, МНМ-оценки обеспечивают приемлемые результаты. Однако при одностороннем характере выбросов наблюдается смещение оценок [10]. Для оценивания моделей (2), (4), (5) в условиях выбросов МНМ уже не обеспечивает состоятельности оценок [8]. В качестве альтернативы в [11] предложен обобщенный метод наименьших модулей (ОМНМ). Без потери общности для удобства будем рассматривать далее модель (1), так как остальные линейные модели формально можно свести к выражению (1).
ОМНМ-оценки для задачи (1) находят как решение задачи
Щ(а) = Е>г (( -хга|) тп, (7)
где р() - некоторая монотонно возрастающая, дважды непрерывно-дифференцируемая на положительной полуоси функция, причем р(0) = 0 ,
vt > о о < р\г) , < р'(0 < о.
МНМ и ОМНМ объединяет общее свойство - в обоих случаях решение находится в узловых точках [9, 12].
Введем гиперплоскости Ц,. = О(а, х,, у) в виде уравнений
у-(а, х) = 0, (. = 1,2,..., п). (8)
Зададим также узловые точки пересечения гиперплоскостей (8):
кт
а(к,,...,кт) = п Ц, 1 < К < к2 <.< кт < п, к, е N. (9)
/=к1
Обозначим как Ц = {Ц1, Ц2,..., Цп} - множество всех гиперплоскостей, и - множество всех узловых точек (9). Поскольку минимум целевых функций Q(a) и Ща) всегда принадлежит множеству и, поэтому точное решение а* задач (6) и (7) может быть получено путем перебора всех узловых точек и выбора в качестве решения той, которая обеспечивает минимум целевой функции.
Перебор всех узловых точек требует решения Ст систем линейных уравнений порядка т и вычисления значения целевой функции в каждой узловой точке а(к[ к ). Это означает, что с ростом п и т наблюдается экспоненциальный рост вычислительных затрат. Фактически практическое применение переборного алгоритма ограничено объемом выборки п < 150 и числом коэффициентов регрессии т < 4 .
В настоящее время предложен ряд других методов решения задач (6) и (7). Однако они также для получения точного или близкого к точному решения требуют больших вычислительных затрат.
Постановка задачи
Целью статьи является исследование особенностей задач (6) и (7) для выяснения причин недостаточно высокой скорости сходимости приближенных методов их решения и проблем, связанных с получением точных реше-
ний, а также описание эффективных вычислительных алгоритмов нахождения МНМ- и ОМНМ-оценок линейных моделей.
Обзор методов нахождения МНМ- и ОМНМ-оценок линейных моделей
Эффективный алгоритм поиска минимума целевой функции должен учитывать ее свойства и особенности геометрии. Исследования показали, что у функции Q(a), несмотря на кусочно-линейный вид, очень «плохая» геометрия -она имеет большое число «оврагов» (равно Ст-1) в виде прямых, причем эти овраги вытянуты, скрещены и многие из них почти параллельны и находятся на близком расстоянии друг от друга, образуя так называемый «пучок». При этом минимум находится внутри этого пучка. Стенки этих оврагов представляют собой выпуклые линейные гиперграни, которые постепенно при приближении к минимуму целевой функции также становятся почти параллельными. В качестве иллюстрации для случая для т = 2, п = 128 и нормально распре деленных случайных ошибок на рис. 1 показаны прямые у{ — а1 — а2х = 0, , = 1,2,..., п, на которых расположены все овраги функции Q(a).
Рис. 1. Графики прямых, на которых расположены все овраги целевой функции Q(a), для случая т = 2 и п = 128 и нормально распределенных случайных ошибок
Указанная специфика геометрии целевой функции приводит ко многим затруднениям при реализации разных алгоритмов нахождения точного решения задачи (6) на основе идей линейного программирования [7, 13-17]. Все эти алгоритмы основаны на симплекс-методе либо на различных его модифи-
кациях. Они осуществляют движение по узловым точкам, являющимся базисными решениями.
Но решение задачи линейного программирования ограничено случаем лишь нескольких десятков наблюдений. Это вызвано накоплением погрешностей из-за ошибок округлений. В результате этого процесс движения по базисным решениям может не дать точного решения задачи (6) или даже зациклиться. Кроме того, требуется хранить в оперативной памяти большой массив данных в виде симплекс-таблиц. Использование техники безошибочных вычислений [18] при устранении роста вычислительных погрешностей приведет к существенному усложнению реализации и росту вычислительных затрат.
При переходе от одной к другой узловой точке указанные алгоритмы не учитывают эту «плохую» геометрию целевой функции.
Объясним, почему переходы по симплекс-таблицам не эффективны. С учетом специфики самой задачи, а именно того, что овраги целевой функции вытянутые, скрещены и многие из них почти параллельны и находятся на близком расстоянии друг от друга, такие переходы без учета геометрии целевой функции часто могут отдаляться от пучка этих гиперплоскостей вида (8), внутри которого лежит точное решение. После этого появляется потребность к возвращению к этому пучку гиперплоскостей, а конкретнее, к точному решению, и для этого осуществляются какие-то спиралевидные неэффективные переходы от одной узловой точки к другой. То есть процесс затягивается, кроме того, при этих спиралевидных переходах также накапливаются вычислительные погрешности. Все эти проблемы существенно усиливаются с увеличением размерности модели и объема выборки.
Известные алгоритмы недостаточно используют свойство ее выпуклости. А именно, ввиду выпуклости целевой функции спуск к точному решению можно осуществлять по различным направлениям. Поэтому нужно выбирать такие направления спуска, которые не позволяли бы отдаляться от пучка гиперплоскостей.
Из приближенных численных методов решения задачи (6) можно отметить метод вариационно-взвешенных квадратических приближений (называемый также алгоритмом Вейсфельда) [7, 19] и методы спуска нулевого порядка [9]. Метод вариационно-взвешенных квадратических приближений не обладает свойством сходимости к точному решению из-за неограниченного возрастания нескольких весовых коэффициентов, а предложенная регуляризация существенно ограничивает точность приближенного решения. Кроме того, с ростом тип резко снижается скорость сходимости алгоритма в окрестности минимума целевой функции. Методы нулевого порядка просты и эффективны с вычислительной точки зрения. Однако они оказались чувствительны к «плохой» геометрии целевой функции, особенно при увеличении размерности модели т. В процессе перехода от одного приближенного решения к другому происходит перескок на другой овраг, который может проходить в стороне от точного решения, что с учетом близости к нулю градиента целевой функции в точках дифференцируемости часто приводит к несостоятельности оценок свободного члена, а значит, и вектора а в целом. Чтобы избежать этого, необходимо значительно повысить точность алгоритма, а это приведет к росту вычислительных затрат.
Анализ известных решений показал, что их основной недостаток состоит в неучете овражности и близости к нулю градиента целевой функции в точках дифференцируемости в достаточно большой окрестности ее минимума. Поскольку целевая функция выпуклая, то спуск у овражной функции наиболее эффективен по ее оврагам. Но этими оврагами являются отрезки прямых (назовем их узловыми прямыми), образованных пересечением гиперплоскостей (8). Таким образом, спускаясь по узловым прямым, мы устраняем все указанные выше проблемы известных методов решения задачи (6).
Целевая функция Ж(а) также имеет овраги, расположенные вдоль узловых прямых. Кроме того, она является вогнутой (выпуклой вверх) и имеет множество локальных минимумов, расположенных в узловых точках, что существенно затрудняет нахождение точного решения задачи (7). Поэтому итерационный алгоритм [20, 21], основанный на идеях линейного программирования, весьма трудоемок и не гарантирует сходимости к точному решению задачи. В [22] для уменьшения вычислительных затрат предложена модификация ОМНМ, основанная на разбиении исходной выборки на подвы-борки и формировании множества оценок по всем подвыборкам. В данном случае повышение точности оценивания достигается за счет осреднения полученных оценок и повторного случайного формирования подвыборок. Это приводит к росту вычислительных затрат. Данный алгоритм также не гарантирует получения точного решения задачи (7).
У целевой функции Ж(а) глобальный минимум также находится в узловой точке, и она имеет вытянутые овраги, расположенные на узловых прямых. Поэтому поиск глобального минимума функции Ж(а) также следует осуществлять по узловым прямым.
Обычный спуск. Данный алгоритм точного решения задачи (6) основан на спуске к точному решению, двигаясь вдоль узловых прямых /^ к ^, каждая из которых является пересечением (т — 1) различных гиперплоскостей Д-:
В качестве начального приближения берется узловая точка а^ к ), являющаяся пересечением т произвольных различных гиперплоскостей Ок.,...,Ок . Исключив одну из гиперплоскостей, получим узловую прямую
/(д. к ). Через любую узловую точку проходит т узловых прямых. Выберем
из этих узловых прямых ту, вдоль которой целевая функция достигает наименьшего значения, которое всегда будет достигаться в одной из узловых точек. Найдя эту узловую точку, продолжим движение из нее по тому же принципу. В результате будет найдена узловая точка а(^ к ) = а*, спуск из
которой невозможен. Эта узловая точка будет являться точным решением задачи (6). В основе алгоритма лежат следующие теоремы.
Реализация МНМ на основе спуска по узловым прямым
Теорема 1. Рассмотрим модель (1), для которой имеется выборка наблюдений (х,, у1) = (хй, хл2,..., хы, у ), (/ = 1, ... , п), и пусть заданы функция Q(a) задачи (6), гиперплоскости (3) и множество и всех узловых точек (9). Тогда функция (2) всегда имеет точку глобального минимума, эта точка либо единственна и принадлежит и, либо состоит из выпуклого линейного многогранника, вершины которого являются точками из и.
Доказательство теоремы 1. Известно, что выпуклая непрерывная кусочно-линейная функция либо имеет глобальный минимум, либо стремится к минус бесконечности. А поскольку функция Q(a) является также ограниченной снизу (Q(a) > 0 , как сумма модулей) функцией, то она всегда имеет точку глобального минимума.
Пусть а* = (а*,а*,...,а*т)Т — стационарная точка и является точкой минимума функции Q(a). Тогда ее градиент в этой точке равен нулю. Поскольку Q(a) - кусочно-линейная функция, то из равенства grad Q(a*) = 0 следует, что функция Q(a) является постоянной функцией на выпуклом многограннике с вершинами а1, а2,..., а1, гранями которого являются гиперплоскости (8), имеющие вид [23]
А = {а: а = £ X ¡ак, £ X = 1}.
к =1 к=1
Функция Q(a) является постоянной до граничных точек многогранника А, лежащих на гиперплоскостях (3). Следовательно, вершины а1, а2,..., а1 этого многогранника являются узловыми точками и функция Q(a) достигает минимума и в этих точках.
Все гиперплоскости (9), и только они, являются особыми точками функции Q(a), поскольку только в них она не дифференцируема. Если взять (т — 1) произвольных невырожденных соотношений вида ^ = у { — (а, х^ = 0, то они в совокупности определяют прямую
кт—1
1(^.1): ПО, к1 е {1,2,...,п} (10)
/=¿1
в пространстве (а1,а2,...,ат) и вместе с тем плоскость Р^ к параллельную оси Q в пространстве а1,а2,...,ат). Присоединяя к системе (10) выражение (6) и рассматривая их совместно, найдем уравнение ломаной М, полученной в результате пересечения поверхности (6) плоскостью Р^ ^ к ).
Если с помощью уравнений, входящих в систему (10), выразить (т — 1) неизвестных а1,а2,...,ат—1 через оставшееся неизвестное и подставить в выражение (6), то получим уравнение проекции ломаной М на плоскость ат) (рис. 2). Точки Т/ ,Т/ ,... ,Т/ являются проекциями на эту плоскость точек пересечения прямой к ) гиперплоскостями (8), не вошедшими в (10).
си
Г
\
К
'1 ^'г
Т,
Рис. 2. Вид функции Q(am)
Функция Q(am) на плоскости (Q, ат) - выпуклая и кусочно-линейная с особыми точками Т,, Т2,...,Т которые являются узловыми. Отсюда следует, что функция Q(a) достигает минимума в узловой точке.
Теорема 2. Алгоритм спуска вдоль узловых прямых для нахождения решения задачи (6) сходится к точному решению за конечное число шагов.
Доказательство теоремы 2. На каждом шаге алгоритма мы находим узловую точку с все меньшим значением целевой функции, а поскольку количество узловых точек конечно, алгоритм выполняется за конечное число шагов.
Докажем, что алгоритм всегда будет достигать минимума целевой функции (т.е. решения). Допустим, что мы находимся на к-м шаге алгоритма, в какой-то узловой точке и пытаемся осуществить спуск вдоль одной из т узловых прямых, проходящих через эту точку, к другой узловой точке с меньшим значением целевой функции. Покажем, что либо мы сможем найти такую точку, либо текущая узловая точка является минимальной.
Целевая функция Q(a) является выпуклой функцией в пространстве Rm. Поэтому, если взять произвольную точку х и т линейно независимых прямых /j (] = 1,2,..., т), проходящих через нее, либо найдется такая окрестность этой точки, в которой целевая функция будет убывать хотя бы вдоль одного из направлений, либо в текущей узловой точке достигается минимум целевой функции.
Поскольку вероятность того, что случайная величина (случайные погрешности) более одного раза примет одно и то же значение равна нулю, то узловые прямые, проходящие через узловую точку, являются линейно независимыми. То есть вдоль одной из них целевая функция убывает, причем она будет убывать вплоть до следующей узловой точки на этой узловой прямой, поскольку между двумя соседними узловыми точками целевая функция (2) будет линейна, так как все ее подмодульные выражения не будут менять знак. Действительно, смена знака означает пересечение одной из образующих гиперплоскостей, что приведет к образованию узловой точки, а мы рассматриваем часть узловой прямой между двумя соседними узловыми точками.
На практике из-за конечной точности измерений гипотетически могут возникнуть ситуации, когда узловые прямые, проходящие через узловую точ-
ку, окажутся параллельными. В этом случае узловой точки не существует, и переходим к рассмотрению другой узловой прямой.
Таким образом, либо произвольная узловая точка является минимумом целевой функции, либо хотя бы по одной из проходящих через нее узловых прямых можно сделать переход к узловой точке с меньшим значением целевой функции. Теорема доказана.
Отметим, что вычислительные затраты данного алгоритма можно значительно сократить, если использовать информацию о предыдущей найденной узловой точке.
Спуск с использованием разреженных матриц. Двигаясь вдоль прямой к 1) для нахождения узловых точек, принадлежащих этой прямой, нужно
для каждой точки решать систему линейных алгебраических уравнений (СЛАУ) порядка т:
а1 + а2Хк1,2 + азХ^,3 + ... + а тХк1 ,т = У^ , а1 + а2 Хк, ,2 + а3 Хк2,3 + ... + атХк, т = Ук,
а. + а2Х, 2 + а3Х, 3 +... + а Х, = у,
1 2 кт—1,2 3 кт—1,3 т кт—1,т ^ кт
а + а2Х.2 + а3Х.3 +... + а Х. = у.,
1 2 .,2 3 .,3 т г,т * .'
(11)
где 1 <, <к2 <•..<кт—1 <п , .е {1,2,.,п}, .€ {к^к2,...,кт—1}.
Очевидно, что СЛАУ двух различных узловых точек, принадлежащих этой прямой одной прямой, отличаются лишь одним (последним) уравнением. Следовательно, вычислительная эффективность алгоритма спуска существенно повысится, если для нахождения узловых точек, которые лежат на прямой к 1), первые (т — 1) строк расширенной матрицы, соответствующей СЛАУ
(11), предварительно преобразуем с помощью элементарных преобразований к ступенчатому виду. Расширенная матрица СЛУ прямой к ) имеет вид
А
№,-,к„ —1)
Г 1 Хк1,2 Хк1,3
1 Хк2,2 Хк2,3
1 Хк3,2 Хк3,3
V 1 Хкт —1,2 Хкт —1,3
, ,т—1
Х.
к3 ,т—1
Хкт —1 ,т—1 Хкт.
Уц
Ук2
Ук3
У к т —1
Применив алгоритм прямого хода метода Гаусса, преобразуем матрицу А(к к ) к ступенчатому виду:
А(
(ки...,кт —1)
Г1 Хк1,2 Хк1,3 ••• ХА1, т—1 Хк1 ,т Уц
0 1 / Х'2,3 ••• / к2, т—1 / Х, к2, т у'2
0 0 1 ... / к3 ,т—1 / х, Л3 ,т у13
V 0 0 0 ... 1 / Хкт —1,т У'т—
197
Л
Используя ступенчатую матрицу А'^ к ^, можно значительно сократить вычислительные затраты на нахождение всех узловых точек, лежащих на прямой ^ к ). Действительно для каждой искомой узловой точки имеем расширенную матрицу
Хк1,т-1 Хк1,т Ук1 / / /
Хк2 ,т-1 Хк2 ,т Ук2 ^ ^ ^
Хк3,т-1 Хк3,т Ук3 (12)
1 Хкт-1 ,т Укт-1 Х , Х у
г,т-1 г,т ^ г у
Варьируя номер г в (12), найдем все узловые точки, лежащие на прямой
1№,...,к„-1) .
Вычислительную эффективность алгоритма спуска можно повысить, если учитывать направление спуска.
Спуск с использованием разреженных матриц и с учетом направления спуска. Используя ступенчатую матрицу А'^ к ) и решив СЛАУ, соответствующую расширенной матрице (12), находим значение т-го коэффициента а^К^) , 1 < к1< к2 <...< кт-1 < п, г е {1,2,..., п), г ¿{^ к2,..., кт_1> для каждой узловой точки. После этого по возрастанию а(тк1,к2,'"кт1,г) упорядочиваем все узловые точки, которые лежат на прямой /^ к ), и выполняем описанный выше алгоритм спуска, но с учетом направления. Если при непосредственном переходе от одной узловой точки к другой значение целевой функции увеличивается, то в этом направлении значение целевой функции будет увеличиваться во всех узловых точках (вытекает из выпуклости целевой функции). Назовем такое направление «плохим». Для осуществления спуска до вычисления значения целевой функции в очередной узловой точке рассматриваем направление спуска. Если оно «плохое», то переходим к следующей точке, не вычисляя в данной узловой точке ни значение целевой функции, ни значения коэффициентов а^,*2,'"кга-1,г), у = 1, 2,.,т -1.
В [24] доказано, что алгоритм спуска по узловым прямым имеет вычислительную сложность Ж = й(т2и2 + т4п 1п п + т2п 1п2 п) . Отметим, что для п > тах(1п2 п; т21п п) вычислительная сложность спуска по узловым прямым Ж = 0(т2 п2).
Нахождение ОМНМ-оценок линейных моделей на основе спуска по узловым прямым
Для реализации идеи неполного перебора узловых точек необходима априорная информация о примерном расположении точного решения а*
А
Г1 Хк1,2 Хк1,3
0 1 / Хк2,3
0 0 1
(К—К-^)
0 0
1 Хг2
задачи (7). Многочисленные вычислительные эксперименты методом статистических испытаний Монте-Карло [25] показали, что вектор а*, как правило, мало отличается от решения задачи (7) методом наименьших модулей, равного
п
а(0) = а^тт£|у. —(а,хЛ|. (13)
г=1
Описанный выше алгоритм реализации МНМ на основе спуска по узловым прямым имеет число операций порядка 0(т2п2). Это пренебрежимо мало
по сравнению с затратами на реализацию полного перебора узловых точек.
Алгоритм решения задачи (7) также основан на использовании узловых прямых и состоит в следующем.
1. Используем в качестве начального приближения а(1) узловую точку решения а(0) задачи (13) методом наименьших модулей, т.е. считаем а(1) = а(0).
2. Уточняем решение путем спуска из узловой точки а(1), поочередно двигаясь вдоль каждой из т узловых прямых, проходящих через а(1). Выберем ту узловую точку, в которой целевая функция Ща) достигает наименьшего значения. Найдя эту точку, продолжим движение из нее по тому же принципу. В результате будет найдена узловая точка а(2), спуск из которой невозможен. При движении вдоль узловой прямой используем ступенчатые расширенные матрицы вида (12).
3. Для найденной узловой точки а(2) определим множество гиперплоскостей О(а(2)) с О, состоящее из а -100% гиперплоскостей, наименее удаленных от а(2) .
Каждая из гиперплоскостей (8) включает в себя случайную ошибку. Ошибки группируются относительно нуля, чем меньше ошибка, тем меньше гиперплоскости отличаются друг от друга. Решение обычно лежит в узловых точках, образованных гиперплоскостями с малыми величинами ошибок. Все эти гиперплоскости будут располагаться в пространстве очень «тесно».
4. С помощью перебора среди всех узловых точек, получаемых в результате пересечения гиперплоскостей из множества О(а(2)), находим ту точку а(3), в которой целевая функция принимает наименьшее значение.
5. Найденную узловую точку а(3) сравниваем с начальным приближением а(1). Если они совпадают, то работа алгоритма останавливается и в качестве решения задачи (2) берется полученная узловая точка. В противном случае в качестве начального приближения а(1) берем полученную узловую точку, переходим к п. 2 и повторяем вышеописанные действия.
Пример. Для нескольких значений п и т проведем исследование описанного алгоритма. Считаем для определенности, что случайная ошибка распределена по стандартному нормальному закону распределения. Воспользуемся методом Монте-Карло. Число испытаний зададим равным 1000. В качестве функции потерь используем р(|х|) = агС§|х| .
В табл. 1 приведено отношение вычислительных затрат переборного алгоритма к затратам предложенного алгоритма спуска по узловым прямым.
Таблица 1
Отношение вычислительных затрат переборного алгоритма к затратам предложенного алгоритма спуска по узловым прямым
n m = 3 m = 4 m = 5
32 7,2 40,4 202,5
64 46,0 414,4 3961,6
128 158,4 1157,9 6203,2
256 222,1 1273,4 6875,0
512 225,4 1339,3 7363,6
Алгоритм спуска по узловым прямым показал высокую точность. В частности, для m = 3, n = 128 в 995 случаях из 1000 полученное решение совпало с точным решением по переборному алгоритму. В остальных пяти случаях средняя относительная погрешность по сравнению с точным решением составила 0,11 %.
Заключение
1. Предложены эффективные алгоритмы реализации метода наименьших модулей и обобщенного метода наименьших модулей при оценивании параметров линейных моделей, основанные на спуске по узловым прямым.
2. Данные алгоритмы позволяют значительно снизить вычислительные затраты по сравнению с переборным алгоритмом.
3. Вычислительная сложность спуска по узловым прямым позволяет на практике реализовать предложенные алгоритмы для анализа экспериментальных данных.
Библиографический список
1. Айвазян, С. А. Прикладная статистика: исследование зависимостей / С. А. Айвазян, И. С. Енюков, Л. Д. Мешалкин. - М. : Финансы и статистика, 1985. - 488 с.
2. Суслов, В. И. Эконометрия / В. И. Суслов, Н. М. Ибрагимов, Л. П. Талышева,
A. А. Цыплаков. - Новосибирск : Изд-во СО РАН, 2005. - 744 с.
3. Кендалл, М. Статистические выводы и связи : пер. с англ. / М. Кендалл, А. Стью-арт. - М. : Наука : ФИЗМАТЛИТ, 1973. - 900 с.
4. Крамер, Г. Математические методы статистики / Г. Крамер ; пер. с англ. - М. : Мир, 1975. - 648 с.
5. Демиденко, Е. З. Линейная и нелинейная регрессия / Е. З. Демиденко. - М. : Финансы и статистика, 1981. - 302 с.
6. Bloomfield, P. Least absolute seviations: theory, applications, and algorithms / P. Bloomfield, W. L. Steiger. - Boston-Basel-Stuttgart : Birkhauser, 1983. - 349 p.
7. Мудров, В. И. Методы обработки измерений. Квазиправдоподобные оценки /
B. И. Мудров, В. Л. Кушко. - М. : Радио и связь, 1983. - 304 с.
8. Болдин, М. В. Знаковый статистический анализ линейных моделей / М. В. Болдин, Г. И. Симонова, Ю. Н. Тюрин. - М. : Наука : ФИЗМАТЛИТ, 1997. - 288 с.
9. Тырсин, А. Н. Оценивание линейных регрессионных уравнений с помощью метода наименьших модулей / А. Н. Тырсин, К. Е. Максимов // Заводская лаборатория. Диагностика материалов. - 2012. - Т. 78, № 7. - С. 65-71.
10. Тырсин, А. Н. Робастное построение линейных регрессионных моделей по экспериментальным данным / А. Н. Тырсин // Заводская лаборатория. Диагностика материалов. - 2005. - Т. 71, № 11. - С. 53-57.
11. Тырсин, А. Н. Робастное построение регрессионных зависимостей на основе обобщенного метода наименьших модулей / А. Н. Тырсин // Записки научных семинаров ПОМИ. - 2005. - Т. 328. - С. 236-250.
12. Тырсин, А. Н. Оценивание линейной регрессии на основе обобщенного метода наименьших модулей / А. Н. Тырсин, Л. А. Соколов // Вестник Самарского государственного технического университета. Серия Физико-математические науки. -2010. - № 5 (21). - С. 134-142.
13. Fisher, W. D. A note on curve fitting with minimum deviations by linear programming / W. D. Fisher // Journal of the American Statistical Association. - 1961. - Vol. 56. -P. 359-362.
14. Barrodale, I. An improved algorithm for discrete L1 linear approximation / I. Barrodale, F. D. K. Roberts // SIAM Journal on Numerical Analysis - 1973. - Vol. 10. -P. 839-848.
15. Narula, S. C. Algorithm AS108: Multiple linear regression with minimum sum of absolute errors / S. C. Narula, J. F. Wellington // Applied Statistics. - 1977. - Vol. 26. -P. 106-111.
16. Armstrong, R. D. Algorithm AS132: Least absolute value estimates for a simple linear regression problem / R. D. Armstrong, D. S. Kung // Applied Statistics. - 1978. -Vol. 27. - P. 363-366.
17. Wesolowsky, G. O. A new descent algorithm for the least absolute value regression problem / G. O. Wesolowsky // Communications in Statistics, Simulation and Computation. - 1981. - Vol. B10, № 5. - P. 479-491.
18. Грегори, Р. Безошибочные вычисления. Методы и приложения / Р. Грегори, Е. Кришнамурти. - М. : Мир, 1988. - 208 с.
19. Weiszfeld, E. On the point for which the sum of the distances to n given points is minimum / E. Weiszfeld // Annals of Operations Research. - 2008. - Vol. 167. -P. 7-41.
20. Панюков, А. В. Взаимосвязь взвешенного и обобщенного вариантов метода наименьших модулей / А. В. Панюков, А. Н. Тырсин // Известия Челябинского научного центра. - 2007. - № 1 (35). - C. 6-11.
21. Панюков, А. В. Об устойчивом оценивании параметров авторегрессионных моделей на основе обобщенного варианта метода наименьших модулей / А. В. Паню-ков // Вестник НГУЭУ. - 2015. - № 4. - C. 339-346.
22. Тырсин, А. Н. Об одном алгоритме реализации обобщенного метода наименьших модулей / А. Н. Тырсин, А. А. Азарян // Обозрение прикладной и промышленной математики. - 2017. - Т. 24, В. 4. - С. 375-376.
23. Рокафеллар, Р. Выпуклый анализ : пер. с англ. / Р. Рокафеллар. - М. : Мир, 1973. -470 с.
24. Тырсин, А. Н. Точные алгоритмы реализации метода наименьших модулей на основе спуска по узловым прямым / А. Н. Тырсин, А. А. Азарян // Вестник Бурятского государственного университета. Математика, информатика. - 2017. - № 4. -С. 21-32.
25. Ермаков, С. М. Метод Монте-Карло и смежные вопросы / С. М. Ермаков. -2-е изд., доп. - М. : Наука : ФИЗМАТЛИТ, 1975. - 472 с.
Тырсин Александр Николаевич
доктор технических наук, доцент, заведующий кафедрой прикладной математики,
Уральский энергетический институт, Уральский федеральный университет им. первого Президента России Б. Н. Ельцина E-mail: [email protected]
Азарян Алексан Артурович аспирант,
Уральский энергетический институт, Уральский федеральный университет им. первого Президента России Б. Н. Ельцина
E-mail: [email protected]
Tyrsin Alexander Nikolayevich doctor of technical sciences, associate professor, head of sub-department of applied mathematics, Ural power engineering Institute, Ural Federal University named after the first President of Russia B. N. Yeltsin
Azaryan Alexan Arturovich postgraduate student, Ural power engineering Institute, Ural Federal University named after the first President of Russia B. N. Yeltsin
УДК 519.237.5:519.24 Тырсин, А. Н.
Методы устойчивого построения линейных моделей на основе спуска по узловым прямым / А. Н. Тырсин, А. А. Азарян // Модели, системы, сети в экономике, технике, природе и обществе. - 2018. - № 1 (25). - С. 188-202.
