function [amax,umax,envs,ccod] = lintregs(A,b,varargin)
// Построение линейной интервальной регрессии.
// LINTREGS(A,b) выдаёт значения коэффициентов линейной регрессионной
// зависимости по матрице входных значений A и интервальному вектору
// выходных данных b путём вычисления максимума распознающего функционала
// множества решений интервальной линейной системы Ax = b.
//
// Синтаксис вызова:
// [amax,umax,envs,ccod] = lintregs(A,b,iprn,epsf,epsx,epsg,maxitn)
//
// Обязательные входные аргументы функции:
// A - точечная mxn-матрица интервальной системы линейных
// алгебраических уравнений, в общем случае прямоугольна;
// b - mx2-вектор нижних и верхних границ правой части
// интервальной системы линейных алгебраических уравнений.
//
// Необязательные входные аргументы функции:
// iprn - выдача протокола поиска; если iprn > 0 - информация
// о ходе процесса печатается через каждые iprn-итераций;
// если iprn <= 0 (значение по умолчанию), печати нет;
// epsf - допуск на точность по значению целевого функционала,
// по умолчанию устанавливается 1.e-6;
// epsx - допуск на точность по аргументу целевого функционала,
// по умолчанию устанавливается 1.e-6;
// epsg - допуск на малость нормы суперградиента функционала,
// по умолчанию устанавливается 1.e-6;
// maxitn - ограничение на количество шагов алгоритма,
// по умолчанию устанавливается 2000.
//
// Выходные аргументы функции:
// umax - значение максимума распознающего функционала;
// amax - доставляющий его вектор значений аргумента,который
// лежит в информационном множестве задачи при umax>=0;
// envs - вектор значений образующих распознающего функционала
// в точке максимума, упорядоченный по возрастанию;
// ccod - код завершения алгоритма (1 - по допуску epsf на
// изменения значений функционала, 2 - по допуску epsg
// на суперградиент, 3 - по допуску epsx на вариацию
// аргумента, 4 - по числу итераций, 5 - не найден
// максимум по направлению).
////////////////////////////////////////////////////////////////////////////////
//
// Эта программа выполняет исследование разрешимости интервальной линейной
// системы уравнений Ax = b с точечной матрицей A и интервальным вектором
// правой части b с помощью максимизации распознающего функционала Uni
// объединённого множества решений этой системы. См. подробности в
// Шарый С.П. Конечномерный интервальный анализ. - Новосибирск: XYZ,
// 2010. - Электронная книга, доступная на http://www.nsc.ru/interval,
// параграф 5.5;
//
// Для максимизации вогнутого распознающего функционала используется
// вариант алгоритма суперградиентного подъёма с растяжением пространства
// в направлении разности последовательных суперградиентов, предложенный
// (для случая минимизации) в работе
// Шор Н.З., Журбенко Н.Г. Метод минимизации, использующий операцию
// растяжения пространства в направлении разности двух последовательных
// градинетов // Кибернетика. - 1971. - №3. - С. 51-59.
// В качестве основы этой части программы использована процедура негладкой
// оптимизации ralgb5, разработанная и реализованная П.И.Стецюком (Институт
// кибернетики НАН Украины, Киев).
//
// С.П. Шарый, 2011 год.
//
////////////////////////////////////////////////////////////////////////////////
//
// проверка корректности входных данных
//
m = size(A,1); n = size(A,2); // n - размерность пространства переменных
k = size(b,1); l = size(b,2);
if ( l ~= 2 )
error('Неправильны размеры вектора правой части!')
end
if ( k ~= m )
error('Размеры матрицы не соответствуют размерам правой части!')
end
for i=1:m;
if ( b(i,1) > b(i,2) )
ErrStr = msprintf('%s %s %s','В компоненте', string(i), ...
'интервал правой части неправилен!');
error(ErrStr)
end
end
////////////////////////////////////////////////////////////////////////////////
//
// задание параметров алгоритма суперградиентного подъёма
//
maxitn = 2000; // ограничение на количество шагов алгоритма
nsims = 30; // допустимое количество одинаковых шагов
epsf = 1.e-6; // допуск на изменение значения функционала
epsx = 1.e-6; // допуск на изменение аргумента функционала
epsg = 1.e-6; // допуск на норму суперградиента функционала
alpha = 2.3; // коэффициент растяжения пространства в алгоритме
hs = 1.; // начальная величина шага одномерного поиска
nh = 3; // число одинаковых шагов одномерного поиска
q1 = 1.0; // q1, q2 - параметры адаптивной регулировки
q2 = 1.1; // шагового множителя
iprn = 0; // печать о ходе процесса через каждые iprn-итераций
// (если iprn < 0, то печать подавляется)
////////////////////////////////////////////////////////////////////////////////
//
// переназначение параметров алгоритма, заданных пользователем
//
nargin = argn(2);
if nargin >= 3
iprn = ceil(varargin(1));
if nargin >= 4
epsf = varargin(2);
if nargin >= 5
epsx = varargin(3);
if nargin >= 6
epsg = varargin(4);
if nargin >= 7
maxitn = varargin(5);
end
end
end
end
end
////////////////////////////////////////////////////////////////////////////////
function [f,g,tt] = calcfg(x)
//
// функция, вычисляющая значения f максимизируемого распознающего
// функционала объединённого множества решений и его суперградиента g;
// кроме того, выводится вектор tt значений образующих функционала
//
////////////////////////////////////////////////////////////////////////////////
//
// предварительное размещение рабочих массивов
tt = zeros(m,1);
dd = zeros(n,m);
////////////////////////////////////////////////////////////////////////////////
// вычисляем значение распознающего функционала и матрицу dd,
// составленную из суперградиентов его образующих
s = 0.5*(b(:,1) + b(:,2)) - A*x;
for i = 1:m;
// вычисление значения i-ой образующей и её суперградиента
tt(i) = 0.5*(b(i,2) - b(i,1)) - abs(s(i));
di = A(i,:);
if s(i) < 0
di = -di;
end
dd(:,i) = di';
end
// выбираем минимальную образующую функционала
// и присваиваем общий суперградиент g
[f,mc] = min(tt);
g = dd(:,mc);
endfunction
////////////////////////////////////////////////////////////////////////////////
//
// формируем начальное приближение x как решение либо псевдорешение
// 'средней' точечной системы, если она не слишком плохо обусловлена,
// иначе берём начальным приближением нулевой вектор
//
sv = svd(A);
minsv = min(sv);
maxsv = max(sv);
if ( minsv~=0 & maxsv/minsv < 1.e+23 )
x = A\(0.5*(b(:,1) + b(:,2)));
else
x = zeros(n,1);
end
////////////////////////////////////////////////////////////////////////////////
//
// Рабочие массивы:
// B - матрица обратного преобразования пространства аргументов
// g - вектор суперградиента максимизируемого функционала
// g1,g2 - используются для хранения вспомогательных векторов
// vf - вектор приращений функционала на последних nsims шагах алгоритма
B = eye(n,n);
vf = %inf*ones(nsims,1);
////////////////////////////////////////////////////////////////////////////////
//
// Присваивание строк оформления протокола работы программы
TitLine = 'Протокол максимизации распознающего функционала Uni';
HorLine = '-----------------------------------------------------------';
TabLine = ' Шаг Uni(x) Uni(xx) ВычФун/шаг ВычФун';
////////////////////////////////////////////////////////////////////////////////
//
// установка начальных параметров, вывод заголовка протокола
//
w = 1./alpha - 1.;
lp = iprn;
[f,g0,tt] = calcfg(x);
ff = f; xx = x;
cal = 1; ncals = 1;
if iprn > 0
printf('\n%100s',TitLine);
printf('\n%60s',HorLine);
printf('\n%77s',TabLine);
printf('\n%60s',HorLine);
printf('\n%5u%17g%17g%9u%9u',0,f,ff,cal,ncals);
end
////////////////////////////////////////////////////////////////////////////////
//
// основной цикл программы:
// itn - счётчик числа итераций
// xx - найденная точка максимума функционала
// ff - значение функционала в точке максимума
// cal - количество вычислений функционала на текущем шаге
// ncals - общее количество вычислений целевого функционала
//
for itn = 1:maxitn;
pf = ff;
// критерий останова по норме суперградиента
if norm(g0) < epsg
ccod = 2;
break
end
// вычисляем суперградиент в преобразованном пространстве
g1 = B' * g0;
g = B * g1/norm(g1);
normg = norm(g);
// одномерный подъём по направлению dx:
// cal - счётчик шагов одномерного поиска,
// deltax - вариация аргумента в процессе поиска
r = 1;
cal = 0;
deltax = 0;
while ( r > 0. & cal <= 500 )
cal = cal + 1;
x = x + hs*g;
deltax = deltax + hs*normg;
[f, g1,tt] = calcfg(x);
if f > ff
ff = f;
xx = x;
end
if cal > nh
hs = hs*q2;
end
r = g'*g1;
end
// если превышено число шагов одномерного подъёма, то выход
if cal >= 500
ccod = 5;
break;
end
// если одномерный подъём занял один шаг, уменьшаем величину шага
if cal == 1
hs = hs*q1;
end
// уточняем статистику и при необходимости выводим её
ncals = ncals + cal;
if itn==lp
printf('\n%5u%17g%17g%9u%9u',itn,f,ff,cal,ncals);
lp = lp + iprn;
end
// если вариация аргумента в одномерном поиске мала, то выход
if deltax < epsx
ccod = 3;
break;
end
// пересчитываем матрицу преобразования пространства
dg = B' * (g1 - g0);
xi = dg / norm(dg);
B = B + w*B*xi*xi';
g0 = g1;
// проверка изменения значения функционала на последних nsims шагах,
// относительного либо абсолютного
vf(2:nsims) = vf(1:nsims-1);
vf(1) = abs(ff - pf);
if abs(ff) > 1
deltaf = sum(vf)/abs(ff);
else
deltaf = sum(vf);
end
if deltaf < epsf
ccod = 1; break
end
ccod = 4;
end
umax = ff;
amax = xx;
// переупорядочение образующих по возрастанию
tt = [(1:m)', tt];
[z,ind] = gsort(tt(:,2),'g','i');
envs = tt(ind,:);
////////////////////////////////////////////////////////////////////////////////
//
// вывод результатов работы
if iprn > 0
if modulo(itn,iprn)~=0
printf('\n%5u%17g%17g%9u%9u',itn,f,ff,cal,ncals);
end
printf('\n%60s\n',HorLine);
end
////////////////////////////////////////////////////////////////////////////////
if umax >= 0
printf('\n%48s\n\n','Информационное множество непусто')
else
printf('\n%54s\n\n','Информационное множество пусто')
end
////////////////////////////////////////////////////////////////////////////////
if ( umax < 0. & abs(umax/epsf) < 10 )
printf('%55s\n','Найденное значение максимума сравнимо с заданной точностью')
printf('%50s\n','- можно перезапустить программу с меньшими epsf и epsx,')
printf('%55s\n\n','чтобы получить больше информации о совместности системы')
end
////////////////////////////////////////////////////////////////////////////////
endfunction
////////////////////////////////////////////////////////////////////////////////