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

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 targ2a_forum
% проба совмещения двух выходов
% сеть: вход [x 1] один слой 4 нейрона
% выход два сумматора Y1 Y2
% веса на входе u, на выходе - v
n = 4; E = eye(4*n);
[x y u v h i] = start(n);
while 1
  [J dE] = grad_jac_n(u,v,x,y);
  [er u v] = step(E,dE,J,u,v,x,y,n);
  i = i + 1;
  h = draw(u, v, x, h, i, er);
  if er < 10^(-5) || 1000 < i
    break
  end
end
end
function [J dE] = grad_jac_n(u,v,x,y)
[J1 J2] = jacobi_mat2(x, u, v);
f = fnet(x, u, v); d = f - y;
J = [J1; J2];
dE = J'*d(:);
end
function [er1 u v] = step(E,dE,J,u,v,x,y,n)
er = sumerr(x, y, u, v);
J2 = J'*J; t = 0; m = 0.0001;
while 1
  H = J2 + m*E; t = t + 1;
  H1 = H^(-1); dw = H1*dE; 
  du = reshape(dw(1:2*n),n,2); 
  dv = reshape(dw(2*n+1:4*n),n,2);
  u1  = u - du'; v1 = v - dv;
  er1 = sumerr(x, y, u1, v1);
  if er < er1
    m = 2*m;
  else
    u = u1; v = v1; break
  end
end
end
function [J1 J2] = jacobi_mat2(x, u, v)
% матрица Якоби для ф-ции:
% y1 = v11*f1 + v21*f2 + v31*f3 + v41*f4
% y2 = v12*f1 + v22*f2 + v32*f3 + v42*f4
n = numel(x); o = ones(n,1); x = [x o];
J1 = [];
for i = 1:size(x,1)
  x1 = x(i,1);
  [du dv] = one_row1(x1, u, v);
  row = [du dv];
  J1 = [J1; row];
end
 
J2 = [];
for i = 1:size(x,1)
  x1 = x(i,1);
  [du dv] = one_row2(x1, u, v);
  row = [du dv];
  J2 = [J2; row];
end
end
function [du dv] = one_row1(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;
 
dv12 = 0; dv22 = 0;
dv32 = 0; dv42 = 0;
 
dv = [dv11 dv21 dv31 dv41 dv12 dv22 dv32 dv42];
end
function [du dv] = one_row2(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,2)*x*df1;
du12 = v(2,2)*x*df2;
du13 = v(3,2)*x*df3;
du14 = v(4,2)*x*df4;
 
du21 = v(1,2)*df1;
du22 = v(2,2)*df2;
du23 = v(3,2)*df3;
du24 = v(4,2)*df4;
 
du = [du11 du12 du13 du14 du21 du22 du23 du24];
 
dv11 = 0; dv21 = 0;
dv31 = 0; dv41 = 0;
 
dv12 = f1; dv22 = f2;
dv32 = f3; dv42 = f4;
 
dv = [dv11 dv21 dv31 dv41 dv12 dv22 dv32 dv42];
end