Запишем решение прямой задачи явно через компоненты векторов измерений и исходных параметров:

i =gi(x1, ..., xK), i = 1,…,N ,

(4.2.1)

где g_i(…) – некоторые функции, в которые входят также и компоненты вектора UX, однако их далее мы договоримся в явных соотношениях не писать. Подставляя (4.2.1) в выражение для невязки (4.1.6) и рассматривая квадрат невязки R² как функцию переменных x_k, k = 1,…,K, для нахождения ее экстремумов получаем систему уравнений , т.е.

, k = 1,…,K

(4.2.1)

В общем случае нелинейных функций g_i непосредственное нахождение решений системы (4.2.2) и анализ их на минимум невязки весьма сложен. Поэтому рассмотрим, для начала, случай линейных функций g_i, который является основой для дальнейшего обобщения его результатов на нелинейные зависимости. Кроме того, задачи нахождения по МНК параметров линейных зависимостей довольно часто возникают на практике, например, именно такие задачи решались в процессе вторичной обработки результатов самолетных измерений потоков солнечного излучения (раздел 3.2).

В случае линейной зависимости запишем (4.2.1) в виде :

.

(4.2.3)

Здесь коэффициенты g_i0, g_ik вовсе не обязаны быть тождественными константами. Они могут быть функциями, и весьма сложными, компонентов вектора UX, важно лишь то, что, поскольку все параметры данного вектора известны и в рамках решения конкретной обратной задачи фиксированы, коэффициенты g_i0, g_ik являются константами в смысле рассматриваемой зависимости измерений от искомых параметров. Подстановка (4.2.3) в (4.2.2) приводит к системе K линейных алгебраических уравнений с K неизвестными

, k = 1,…,K.

(4.2.4)

Используя определенные выше вектора X (x_k), Y (y_i), вводя вектор G₀ (g_i0) и матрицу G (g_ik), i = 1,…,N, k = 1,…,K перепишем (4.2.4) в матричной форме:

где символом “+” обозначено транспонирование матрицы. При записи (4.2.5) векторы считаются столбцами, а первый индекс матрицы – индексом строки, такого порядка договоримся придерживаться и далее. Умножая обе части (4.2.5) слева на (G⁺G)^-1 получаем искомое решение

Заметим, что матрица (G⁺G) системы уравнений (4.2.5) симметрична и положительно определена (по критерию Сильвестра [11]). Отсюда следует, что решение (4.2.6) существует, единственно (поскольку определитель положительно определенной матрицы больше нуля) и соответствует минимуму невязки (поскольку положительно определенная матрица – ее вторая производная). Часто выражение (4.2.6) называют решением системы линейных уравнений G₀ + GX = Y по МНК, этой терминологией мы пользовались в разделе 3.2 и договоримся пользоваться в дальнейшем. Для квадратной матрицы G (4.2.6) переходит в “обычное” решение X = G^-1(Y-G₀).

При практических вычислениях по формулам (4.2.6) на компьютере для уменьшения возможной погрешности, связанной с накоплением машинной ошибки округлений, здесь и далее рекомендуется следующий стандартный прием нормировки [15]. Обозначим для краткости систему (4.2.5) в виде AX = B и ведем вектор , k = 1,…,K. Переходим к системе A'X' = B', где a'_jk = a_jk / (d_j d_k), b'_k = b_k / d_k, после ее решения X'=(A')^-1B, получаем окончательный результат как x_k = x'_k/d_k. В силу симметричности и положительной определенности матрицы системы A' для ее обращения следует использовать очень эффективный метод квадратного корня [16]. Вычисление сомножителей в (4.2.6) следует производить справа налево, тогда все операции сведутся лишь к умножению вектора на матрицу.

Выше мы считали, что вклад в невязку всех квадратов разностей результатов измерений и вычислений одинаков. Однако, часто желательно учесть индивидуальные особенности этих вкладов. В этом случае используется обобщение МНК – МНК “с весами” [16]. Запишем формулу для невязки (4.1.6) в виде :

(4.2.7)

где w_i > 0 – некий “вес”, приписываемый точке i. Тогда, для линейной зависимости (4.2.3) система уравнений (4.2.4) переходит в

, k = 1,…,K.

(4.2.8)

Для записи системы уравнений (4.2.8) и ее решения в матричной форме придется ввести не вектор, а диагональную матрицу весов W(w_ij), w_ii = w_i, w_i,_ji=0, i = 1,…,N, j = 1,…,N. Тогда матрица системы (4.2.8) запишется как (G+WG), свободный член как G⁺W(Y-G₀), и решение как :

Следует заметить, что при практических вычислениях матрицы и свободного члена удобнее пользоваться явными выражениями из (4.2.8). Смысл введения матрицы весов W станет понятным в следующем разделе. Здесь отметим, что решение задачи МНК не зависит от абсолютных значений весов, то есть умножение всех wi на константу не изменяет значения искомых параметров X, в частности, если все wi равны, то решение (4.2.9) совпадает со случаем решения “без весов” (4.2.6).

В принципе веса wi могут выбираться из различных соображений. Стандартной является ситуация, когда в качестве весов выбирается обратный квадрат среднеквадратической погрешности измерений, т.е. w_i = 1 / s_i², где s_i – СКО измерения y_i. Теоретическое обоснование такого выбора будет дано в следующем разделе. Пока отметим его очевидный смысл: чем больше погрешность измерения, тем и меньше его вклад в невязку и тем меньше требования к близости соответствующих значений y_i и _i. Другим важным частным случаем использования весов является переход к относительной величине невязки, т.е. вычисление в R суммы квадратов не абсолютных, а относительных отклонений y_i от _i. Очевидно, что в этом случае выполняется равенство w_i = 1/y_i². Если вычисляются и значения относительной невязки, и задано относительное СКО точек ряда _i, то получаем т.е. случай минимизации относительной невязки с заданием относительного СКО полностью эквивалентен случаю задания для каждой точки абсолютного СКО при минимизации абсолютной невязки. Отметим, что указанную схему “с весами” мы использовали для учета погрешностей измерений в разделе 3.2.

Искомые в МНК параметры x_k должны быть линейно независимы, иначе матрица системы уравнений (4.2.5) будет вырожденной, и обратная ей матрица существовать не будет. Встречаются, однако, случаи, когда требуется учитывать линейные связи между искомыми параметрами, именно такая ситуация была описана в разделе 3.2. при вторичной обработке результатов зондировок. Запишем в общем виде указанные связи как :

j = 1,…,J.

(4.2.10)

Очевидно, что условия (4.2.10) должны быть линейно независимыми и J < K (иначе следует просто исключить из (4.2.10) линейно зависимые строки, уменьшив значение J). По теореме о базисном миноре [11] в условиях (4.2.10) J линейно независимых столбцов, будем для определенности считать их первыми слева (иначе следует просто перенумеровать компоненты x_k). Разобьем вектор X на две части: X^(J) (x_k), k = 1,…,J и : X^(K-J) (x_k), k =J+1,…,K. Тогда в матричной форме условия (4.2.10) запишутся в виде :

где C₀ (c_j0), C^(J) (c_jk), k = 1,…,J, C^(K-J) (c_jk), k = J + 1,…,K, j = 1,…,J.

Матрица C^(J) является невырожденной, следовательно система (4.2.11) разрешима относительно X^(J)

На основе соотношения (4.2.12) получаем выражение всего вектора X через его независимую часть X^(K-J)

где вектор B₀ b_k0 и матрица B b_kl имеют следующую структуру: b_k0 = ((C^(J))^-1(-C₀))_k и b_kl=((C^(J))^-1(-C^(K-J)))_kl, если k = 1,…,J, l = 1,…,J; b_k0 = 0, b_kk = 1 и b_klk = 0, если k = J + 1,…,K; l=J + 1,…,K. Подставляя (4.2.13) в исходную систему уравнений GX = Y, записывая ее решение по МНК для независимых переменных X^(K-J) и переходя по (4.2.13) вновь ко всему вектору X, окончательно, с учетом общего свойства матриц (GB)⁺ = B⁺G⁺ (и добавлением весов), получаем

Именно с помощью соотношения (4.2.14) решалась в разделе 3.2. система уравнений (3.2.6) при учете ограничений на параметры (3.2.7). Заметим, что в (4.2.14) подлежащая обращению матрица, как, впрочем, и все подобные матрицы, которые встретятся далее, по-прежнему является симметричной и положительно определенной. Для эффективного вычисления произведений нескольких матриц следует использовать описанный выше прием умножения матрицы справа на вектор, выбирая последовательно в качестве таких векторов столбцы последней матрицы произведения.

Перейдем теперь к общему случаю нелинейной зависимости измерений Y от параметров X (4.2.1). Возьмем некоторые начальные значения параметров X₀ (x₀,k) и разложим зависимость (4.2.1) в ряд Тейлора, ограничиваясь только первым (линейным) членом

(4.2.15)

Разность есть (в рассматриваемом приближении) линейная функция разности параметров x_k – x_0k, что позволяет, используя полученное выше решение для случая линейной зависимости, легко построить итерационный алгоритм для нелинейной. Этот стандартный прием сведения нелинейных задач к линейным известен как линеаризация. В матричных обозначениях из (4.2.15) имеем

(X) = G0(X - X0) + (X0),

(4.2.16)

где G₀ – матрица частных производных , i = 1,…,N, k = 1,…,K, вычисленных в точке X₀. Это обозначение матрицы производных удобно тем, что в линейном случае (4.2.3) матрица G, очевидно, имеет точно такой же смысл, следовательно, сохраняется преемственность обозначений. Оператор решения прямой задачи также сохраняет прежнее обозначение G(X,UX), но для краткости договоримся далее писать просто G(X). Применяя итерационно к (4.2.16) рассмотренное выше решение по МНК (4.2.9) получаем ,

где n = 0,1,2 – номер итерации.

При практическом применении формул (4.2.17) возникают три сложности: неопределенность в выборе нулевого приближения, необходимость выработки критерия прерывания итераций и возможный сильный разброс значений искомых величин при последовательных итерациях.

Выбор нулевого приближения X₀ в конкретных задачах обычно осуществляется из физических соображений. Он связан с некоторым “угадыванием решения”, действительно, чем ближе будет нулевое приближение к окончательному решению, тем меньше потребуется итераций и тем лучше будет их сходимость. Обычно в качестве X₀ берут некие априорные средние значения, для атмосферно - оптических задач это обычно средние климатические и т.п. данные. Иногда имеется возможность каким-либо способом получить приближенное решение задачи, пусть даже в грубом приближении, тогда его следует брать в качестве X₀, это обычно существенно повышает эффективность итерационного процесса (4.2.17). Заметим, что в силу нелинейности задачи, решение ее может быть не единственным, т.е. зависеть от конкретного выбора X₀. Более подробно вопросы, связанные с этой зависимостью, мы обсудим в разделе 4.4.

Теоретически в качестве критерия прерывания итераций используется стандартное условие (X_n+1,X_n) < , где (…) – некоторая метрика, а – параметр, имеющий смысл точности решения. Обычно в качестве (…) используется евклидова метрика (4.1.5), поскольку она согласована с метрикой измерений, но возможны и любые другие варианты, например жесткое условие вида  [15]. При практических вычислениях, однако, все оказывается сложнее. Накопление машинной погрешности вычислений, а также возможные особенности поведения самой невязки вблизи точки минимума, приводят к тому, что величина (X_n+1,X_n) перестает уменьшаться с ростом n, следовательно, условие останова итераций не может быть выполнено при слишком малом . Поэтому, чтобы решение не зависело от конкретного выбора e, нередко используют и другие критерии останова. Так эффективным приемом является анализ величины (X_n+1,X_n) как функции n и прекращение итераций, когда ее устойчивое уменьшение сменяется на колебания вокруг некоторого значения. В простейшем варианте в этом случае решение об останове принимается в интерактивном режиме [17]. Другой простой прием – выбор решения, соответствующего минимуму невязки, при фиксированном числе итераций. Заметим, что особенности сходимости итераций в любом случае обусловливаются условиями конкретной задачи и обычно требуют специального изучения в рамках предварительных численных экспериментов [17].

Причину иногда возникающего сильного разброса значений (X_n+1,X_n), т.е. сильной разницы в искомых значениях на двух соседних итерациях, легко понять из физического смысла. Действительно, матрица частных производных Gn зависит от текущего значения X_n и вполне возможно, что X_n попадет в область, где измеряемые величины будут очень слабо зависеть от какого-либо компонента x_n,j. Но это значит, что, не влияя существенно на значения измерений, можно очень сильно варьировать значение x_n,j. Именно так обычно и “поступает” алгоритм (4.2.17). Подходы к устранению этого затруднения весьма разнообразны. Все они основаны на парадоксальной идее “замедления сходимости”, в основе которой лежат приемы, не позволяющие значениям X_n на соседних итерациях различаться слишком сильно (по принципу “тише едешь – дальше будешь”). Мы будем еще не раз возвращаться к этой теме, а пока рассмотрим одну из простейших возможностей.

Запишем исходное уравнение (4.2.16) для итерации номер n и прибавим к обеим его частям произведение G_n(X_n – X₀). Проведя элементарные преобразования и решив полученное уравнение по МНК, вместо (4.2.17) получаем согласно результатам работы [18]:

В алгоритме (4.2.18) “отсчет” всех итераций ведется от X0, что является некоторым препятствием для слишком сильного разброса их значений. По крайней мере, примеры практического использования (4.2.18) [17,18] свидетельствуют о более высокой эффективности, чем (4.2.17), несмотря на больший объем требуемых вычислений. Договоримся в дальнейшем в качестве решения задачи МНК использовать именно алгоритм (4.2.18).

Нередко, исходя из физического смысла, на искомые параметры x₁,…,x_k накладываются дополнительные условия (связи и ограничения), т.е. с математической точки зрения возникает задача поиска не безусловного, а условного экстремума. Эта задача более сложная, а общие методы ее решения, например классический метод неопределенных множителей Лагранжа [13], не всегда вписываются в идеологию МНК.

В некоторых частных случаях удается учесть связи и ограничения на значения искомых параметров, используя специальные приемы. Например, в рассмотренном выше случае линейных связей между параметрами x₁,…,x_k, выражаемых через вектор B₀ и матрицу B, элементарно получаем

Алгоритм (4.1.19) является обобщением алгоритма (4.1.14) на случай нелинейных задач.

Особые трудности возникают, если ограничения на возможные значения параметров записываются в виде неравенств. Например, по своему физическому смыслу практически все параметры (концентрации газов и аэрозольных частиц, альбедо поверхности и т.п.) должны быть неотрицательны. В этом случае достаточно очевидный прием “снятия” ограничений состоит в переходе от самих значений x1,…,xk к их логарифмам [19,20]. Однако, строго говоря, в этом случае, из-за нелинейности логарифма, значения логарифмов параметров, при которых достигается минимум невязки, не обязаны соответствовать значениям самих параметров, обладающих аналогичным свойством. То есть логарифмирование фактически вносит дополнительную погрешность в полученное по МНК решение. Поэтому, несмотря на простоту и заманчивость этого приема, к его применению следует относиться осторожно, исследуя его “плюсы и минусы” применительно к условиям конкретных задач.

Вместе с тем, в настоящее время хорошо известен общий метод, позволяющий приближенно учитывать любые сложные связи и ограничения на восстанавливаемые параметры, – метод штрафных функций [13].

Пусть на искомые параметры наложены J условий связей, которые, не ограничивая общности, можно записать в виде

Функции c_j предполагаются непрерывными и дифференцируемыми во всей области значений аргументов. Заметим, что в отличие от линейного случая (4.2.10), в методе штрафных функций в условиях (4.2.20) соотношения между числом связей J и числом параметров K могут быть любыми, в частности, допустимо J K (и, разумеется, конкретные c_j не обязаны зависеть от всех аргументов сразу). Метод штрафных функций состоит в замене поиска минимума невязки R (4.2.7) поиском минимума следующей величины

(4.2.21)

где h_j – некоторые константы. Идея метода элементарна. Действительно, при точном выполнении условий связей (4.2.20) дополнительная сумма R²_H в (4.2.21) с функциями c_j не вносит никакого вклада в величину собственно невязки R², чем же хуже выполняются эти условия (то есть, чем дальше значения c_j от требуемого нуля), тем больше и вклад дополнительной суммы в общую величину R²_C. Этот вклад является как бы “штрафом” за нарушение условий связи, отсюда и название метода (штрафными функциями называются выражения h_jc_j(x₁,…,x_K)). При поиске минимума R²_C решение будет стремиться к значениям параметров, при которых этот дополнительный вклад условий минимален, то есть к наиболее точному выполнению условий связи (4.2.20). Выбор констант h_j, j = 1,…,J в (4.2.21) в принципе достаточно произволен. Понятно, что чем они больше, тем точнее будут выполнены для решения условия связей (4.2.20). Теоретически h_j должны стремиться к бесконечности [13], но, с практической точки зрения, чем больше h_j, тем более нелинейной будет задача и тем сложнее приспособить к ней вычислительный алгоритм. Поэтому на практике при выборе h_j следует соблюдать определенную меру. Почти всегда все h_j выбирают одинаковыми, т.е. управляют алгоритмом при помощи одного параметра штрафных функций h = h₁ = …=h_j.

Для решения задачи поиска минимума величины (4.2.21) применим линеаризацию: сначала получим решение для линейных функций g_i и линейных же c_j, а потом сведем к нему нелинейный случай. Из системы уравнений вместо (4.2.8) с учетом линейных зависимостей (4.2.3) и (4.2.10) получаем

, k = 1,…,K.

(4.2.22)

Переходя к матричным обозначениям и вводя вектор и матрицу связей C₀ (c_j0), C (c_jk), j=1,…,J, k = 1,…,K, а также аналогичную матрице W диагональную матрицу H (h_j1), h_jj = h²_j, h_j,lj=0, получаем решение системы (4.2.22)

Отметим, что на практике, особенно при равенстве всех hj для вычислений, конечно, следует использовать явные выражения (4.2.22). Для нелинейного случая разложим, как и выше, функции g_i и c_j в ряд Тейлора и ограничимся линейным членом в нем. Получим уравнение для итераций

где G_n – введенная выше матрица частных производных , C_n –аналогичная матрица частных производных , j = 1,…,J, и C(X_n) –вектор значений функций c_j (4.2.20). Все векторы и матрицы вычисляются для аргумента X_n. Применяя к (4.2.24) описанный выше прием улучшения сходимости итераций, а именно, прибавляя к обеим частям (G_n⁺WG_n + C_n⁺HC_n)(X_n -  X₀), окончательно получаем итерационный алгоритм МНК с учетом условий (4.2.20) по методу штрафных функций

Важным моментом использования общего выражения (4.2.25) является то, что в рамках текущей итерации значения параметров на предыдущих итерациях уже определены, следовательно, они могут играть роль констант в штрафных функциях. Потребуем, например, чтобы искомые параметры на текущей итерации не слишком отличались от их значений на предыдущей итерации, т.е. используем обсуждавшийся выше прием “замедления сходимости”. Этому, очевидно, соответствуют условия связи

причем X_n здесь уже не переменная, а константа. В этом случае матрица совпадает с матрицей H, поскольку C_n – единичная матрица. Также по условиям (4.2.26) всегда C(X_n) = 0 и (4.2.25) переходит в предложенный в работе [21] алгоритм с улучшенной сходимостью

В алгоритме (4.2.27) чем больше значения весов h_j, тем ближе значения предыдущей и последующей итераций, таким образом, надлежащим подбором h_j можно добиваться плавной, без “разброса” сходимости итераций.

Упомянем еще об одном частном случае применения метода штрафных функций [22], который может быть востребован при решении обратных задач атмосферной оптики. Нередко встречающимся на практике случаем является ситуация, когда часть искомых параметров (или даже все) могут принимать лишь целочисленные значения. Например, можно рассматривать задачу учета некоторого влияющего на перенос излучения фактора, что можно описать введением в вектор X некого параметра, который должен быть равен только нулю или единице, т.е. “включать” или “выключать” указанный фактор. При малом числе подобных параметров такие задачи можно решать методом перебора всех возможных вариантов их значений, однако, при большом их количестве это становится нереальным. Тогда целесообразно использовать метод штрафных функций, поскольку уравнения связей в этом случае имеют вид x_k(x_k - 1) = 0.

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

j(x1,...,xk)0, j = 1,…,J.

(4.2.28)

Но очевидно, что эти условия (4.2.28) можно приближенно свести к рассмотренным выше условиям связи (4.2.20), если вместо y_j использовать дифференцируемые функции, ведущие себя аналогично, т.е. близкие к нулю, там, где неравенство выполняется и принимающие достаточно большие значения в области, где неравенство не выполняется. В качестве наиболее простых условий связи, обладающих требуемыми свойствами, можно использовать, например, условие:

cj = exp(-hjj(x1,...,xk)),     j = 1,…,J .

(4.2.29)

Константы здесь имеют смысл рассмотренных выше, в соотношениях (4.2.20), поскольку матрица частных производных от (4.2.29) равна

 Поэтому формально для согласования с полученными выше формулами следует всюду положить:

Выбор h_j также достаточно произволен, причем, чем они больше, тем точнее выполняются неравенства, но тем сильнее проявляется нелинейность экспоненты. Очевидно, что метод штрафных функций элементарно переносится и на случай, когда на параметры наложены и связи (4.2.20), и условия (4.2.28).