Проблема с обучением простой нейронной сети методом Левенберга-Марквардта

function leven_parabola1
E = eye(12); % единичная матрица
[x y u v h] = start;
i = 0;
while 1
  [dE J] = grad_jac(u, v, x, y);
  [er u v] = one_step(E, dE, J, u, v, x, y);
  i = i + 1;
  h = draw(u, v, x, h, i, er); 
  if er < 10^(-4) || 50 < i
  % выход из цикла когда или ошибка меньше 0.0001 или число итераций больше 50
    break
  end
end
end

function [x y u v h] = start
x = (0.1:0.05:0.9)'; % столбец входных х
y = (2*x - 1).^2 + 0.1; % парабола, которую надо получить на выходе
plot(x, y, '.'); % выводим параболу на график
u = rand(2,4); v = rand(4,1); % генерация двух матриц весов
y1 = fnet(x, u, v);  % функция сети получает значения у через веса сети
h = plot(x,y1,'-w');  % выводим эти значения на график
% h - указатель на линию x,y1 чтобы потом ее можно было стереть
end
function f = fnet(x, u, v)
% функция сети f = v1*f(xu11+u21)+v2*f(xu12+u22)+v3*f(xu13+u23)+v4*f(xu14+u24);
rows = size(x,1);
for i = 1:rows
  x1 = x(i,1);
  s1 = x1*u(1,1) + u(2,1);
  s2 = x1*u(1,2) + u(2,2);
  s3 = x1*u(1,3) + u(2,3);
  s4 = x1*u(1,4) + u(2,4);
  f1 = smd(s1); f2 = smd(s2);
  f3 = smd(s3); f4 = smd(s4);
  f(i,1) = f1*v(1)+f2*v(2)+f3*v(3)+f4*v(4);
end
end
function f = smd(x)
f = 1./(1+exp(-x)); % сигмоид
end
function [dE J] = grad_jac(u, v, x, y)
J = jacobi_mat(u, v, x); % J - матрица Якоби
y1 = fnet(x, u, v); d = y1 - y; % d - столбец разностей
dE = J'*d; % градиент (J' - символ ' - транспонирование матрицы J)
end

function J = jacobi_mat(u, v, x)
% матрица Якоби для функции: y = v11*f1 + v21*f2 + v31*f3 + v41*f4
J = []; n = numel(x); % n - число элементов в массиве х (в данном случае n = 17)
for i = 1:n
  row = one_row(x(i), u, v);
  J = [J; row];
end
end
function row = one_row(x, u, v)
s1 = x*u(1,1)+u(2,1); [f1 df1] = sgm(s1);
s2 = x*u(1,2)+u(2,2); [f2 df2] = sgm(s2);
s3 = x*u(1,3)+u(2,3); [f3 df3] = sgm(s3);
s4 = x*u(1,4)+u(2,4); [f4 df4] = sgm(s4);
 
du11 = v(1,1)*x*df1;
du12 = v(2,1)*x*df2;
du13 = v(3,1)*x*df3;
du14 = v(4,1)*x*df4;
 
du21 = v(1,1)*df1;
du22 = v(2,1)*df2;
du23 = v(3,1)*df3;
du24 = v(4,1)*df4;
 
du = [du11 du12 du13 du14 du21 du22 du23 du24];
 
dv11 = f1; dv21 = f2; dv31 = f3; dv41 = f4;
dv = [dv11 dv21 dv31 dv41];
row = [du dv];
end
function [f df] = sgm(x)
f = 1./(1+exp(-x)); df = f.*(1-f); % сигмоид и его производная
end

function [er1 u v] = one_step(E, dE, J, u, v, x, y)
er = sumerr(x, y, u, v);
J2 = J'*J; % произведение транспонированной матрицы на себя - аппроксимация матрицы Гессе
m = 0.0001;
while 1
  H = J2 + m*E;
  H1 = H^(-1); % Н1 - обратная матрица к Н
  dw = H1*dE;  % dw - столбец изменения весов
  du = reshape(dw(1:8),4,2); % из первых восьми чисел dw создаем матрицу размером 4х2
  dv = reshape(dw(9:12),4,1); % из элементов с 9 по 12 столбца dw создаем матрицу 4х1
  % обновляем веса
  u1  = u - du'; % здесь транспонирование du необходимо ввиду того что du сейчас 4х2, а u 2х4
  v1 = v - dv;
  er1 = sumerr(x, y, u1, v1); % er1 - новая ошибка при новых весах
  if er < er1
    m = 2*m;
  else
    u = u1; v = v1; break
  end
end
end
function err = sumerr(x, y, u, v)
% err - суммарная ошибка сети при данных u, v
y1 = fnet(x, u, v);
d = y1 - y; err = sum(d.*d);
end

function [x y u v h] = start
x = (0.1:0.05:0.9)'; % входной столбец значений х тот же самый
y1 = (2*x - 1).^2 + 0.1;  % парабола на верхнем выходе
y2 = -3*(x - 0.6).^2 + 0.9; % парабола на нижнем выходе
y = [y1 y2]; % выходной массив у состоит уже из двух столбцов
pot(x,y1,'.'); plot(x,y2,'.k'); % выводим два графика: первый синие точки, второй - черные
 
u = rand(2,4); v = rand(4,2); % матрица весов v состоит уже из двух столбцов
 
f = fnet(x, u, v); % функция сети выдает два столбца
h(1) = plot(x,f(:,1),'-c');
h(2) = plot(x,f(:,2),'-y');
end
function f = fnet(x, u, v)
% y1 = v11*f1 + v21*f2 + v31*f3 + v41*f4
% y2 = v12*f1 + v22*f2 + v32*f3 + v42*f4
rows = size(x,1);
for i = 1:rows
  x1 = x(i,1);
  s1 = x1*u(1,1) + u(2,1);
  s2 = x1*u(1,2) + u(2,2);
  s3 = x1*u(1,3) + u(2,3);
  s4 = x1*u(1,4) + u(2,4);
  f1 = smd(s1);
  f2 = smd(s2);
  f3 = smd(s3);
  f4 = smd(s4);
  f(i,1) = f1*v(1,1)+f2*v(2,1)+f3*v(3,1)+f4*v(4,1);
  f(i,2) = f1*v(1,2)+f2*v(2,2)+f3*v(3,2)+f4*v(4,2);
end
end

function [u v] = update_u_v(u, v)
u = reduce_all(u, 7, 2);
v = reduce_all(v, 7, 2);
end
function u = reduce_all(u, lim, div)
global mt; w = abs(u);
while 1
% ф-ция maxmat - находит максимум (val) матрицы (w)
% и возвращает строку (r) и столбец (с) где этот максимум
  [r c val] = mt.maxmat(w);
  if lim < val
    w(r,c) = 0;
    u(r,c) = u(r,c)/div;
  else
    return
  end
end
end