Lipschitz Delay Estimation for NARX Networks (implementation not working?)

Hi,

I am having a hard time trying to recreate the results of example 3 of the paper "A New Method for Identifying Orders of Input-Output Models for Nonlinear Dynamic Systems", https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4793346&tag=1.

The basic idea is to compute the Lipschiz-Number (based on lipschitz quotient) while successively using more and more regressors (Input/Output-delays). If too less regressors are used, the lipschitz number is very high; if too much regressors are used, the lipschitz number won't change that much compared to the lipschitz number computed for the optimal number of input/output delays.

My implementation:

function qn = LipshitV2(y,U,omax)

% Usage: This function can be used to estimate regressor output and input delays based on 
% Lipschitz-Number
 
% Inputs: 
%   y  is a  horizontal vector of output time series
%   U is matrix with each row representing an input time series
%   omax: maximum timedelay to consider
 
% Output:
% q: vector containing Lipschitz-Numbers, i.e. q(1) = q^1 = q^(1,0,...,0);
 

% Keywords: Embedding Dimension, MIMO model order estimation, non-linear system
% identification
 
% Ref:
% -Wang, Hsu et al. (2009). Minimal model dimension/order determination 
% algorithms for recurrent neural networks
%
% -He, Asada (1993). A New Method for Identifying Orders of Input-Output Models
% for Nonlinear Dynamic Systems

%%

[OUT,N]=size(y);
[IN,N]=size(U);
g = floor(0.01*N);
% Compute the regressor-matrix

RegM = regressorMatrix([y', U'], 1,omax);   % Matrix of Regressors
RegM = RegM(1+omax:end,:);                  % Remove NaNs 
[rRegM, cRegM]=size(RegM);

y=y(1+omax:end);                            % shift y to fit regression matrix 
for n=1:cRegM
qnij=[];                                   % to store all left sides of eq 17
   D=dist(RegM(:,1:n)');                   % euclidean distances
   
   for i=1:rRegM   
       d11 = D(i,1:i-1);                   % left distances
       d12 = D(i,i+1:end);                 % right distances
       d   = [d11,d12];                    % distance vector excluding distance to self
       yj   = [y(1,1:i-1),y(1,i+1:end)];   % y(j|=i)
       qnij(i,:) = abs(y(1,i)-yj)./d;      % eq 12 left    
   end
   %qnr = maxk(qnij(:),g);
   qnr = sort(qnij(:),'descend');
   qnr = qnr(1:g);
   qn(n) = (prod(sqrt(n)*qnr))^(1/g);
end

semilogy(qn)
grid on
ylim([0 1e3])
title('Lipschitz-Plot')
xlabel('order n = l1...lk+p');
ylabel('Lipschitz-Number q^n');

Function to create regressor matrix:

function R = regressorMatrix(IN, tdmin,tdmax)
% Usage:  Create a matrix of lagged (time-shifted) series, 
%         similar to the lagmatrix() function of the econometrics toolbox. 
%         Only positive lags (delays) supported.
%
 
% Inputs: 
%   IN:    A (N by NoR) Matrix of NoR time series with N samples
%   tdmin: Minimum delay (shift of the time series)
%   tdmax: Maximum delay
 
% Output:
%   R:      N by NoR*(tdmax-tdmin+1) Matrix, each column beeing shifted 
%           from tdminto tdmax and NaN-padded.


%%
if (tdmin <0)
    error('Only positive lags (delays) are supported');
end
if (tdmin > tdmax)
    error('Minimum delay is bigger than maximum delay')
end

[N,NoR] = size(IN);

R = NaN(N, NoR*(tdmax+(1-tdmin)));
for k = 1:tdmax+(1-tdmin)
    R(k+(tdmin):end, ((k-1)*NoR+1):(k*NoR)) = IN(1:(end-k-tdmin+1),1:NoR);
end
end

Test-script for the system given in example 3 of the original paper on Lipschitz-Number:

% Test Example 3 of original Paper
N = 2000;
y = zeros(1,N+1);
u = y;
for t = 4:N+1    
    u(t) = (rand-0.5)*2;
    x1 = y(t-1);
    x2 = y(t-2);
    x3 = y(t-3); 
    x4 = u(t-1);
    x5 = u(t-2);
    y(t) = (x1*x2*x3*x5*(x3-1)+x4)/(1 + x1^3 + x2^2);
end

%%
LipshitV2(y,u,5);

Plot:

Plot provided by the paper (Pay attention to the differently scaled y-axis. The values are difficult to read, the y-axis labeling starts at 10^0.) :

The computed Lipschitz-Numbers differ significantly from the numbers stated in the paper. In the Plot, there is no clear "saturated" region, as there is in the paper.

I also tried the function lipschit() of the nnsysid20 matlab toolbox: https://de.mathworks.com/matlabcentral/fileexchange/87-nnsysid. The results differ from both my results and the results stated in the paper. Also here no clear saturated region is visible.

I tried different values for g and normalizing the time series with normalize() and mapminmax(), basicly getting the same results.

I am very grateful for any help!

Cheers,

Kris

Answer 1

[Kris Schweikert]

Solved. In case anybody is facing the same issue: There is a typo in eq. 16 in the original paper. It should be x_1^2 instead of x1_^3.

Answer 2

[YINGYING ZHU]

Thank you very much for the code, it was well worth learning