Multivariate Loewner Framework (mLF)

Overview

The Multivariate Loewner Framework is introduced by A.C. Antoulas, I-V. Gosea and C. Poussot-Vassal in "On the Loewner framework, the Kolmogorov superposition theorem, and the curse of dimensionality" (or arXiv version), in SIAM Review (Research Spotlight), Vol. 67(4), pp. 737-770, November 2025. It allows constructing a $n$-variate rational function approximating a $n$-dimensional tensor (either real or complex valued). It is suited to approximate, from any $n$-dimensional tensor

• $n$-variable static functions
• $n$-variable (parametric) dynamical systems

Contributions claim

• We propose a generalized realization form for rational functions in $n$-variables (for any $n$), which are described in the Lagrange basis;
• We show that the $n$-dimensional Loewner matrix can be written as the solution of a series of cascaded Sylvester equations;
• We demonstrate that the required variables to be determined, i.e. the barycentric coefficients, can be computed using a sequence of small-scale 1-dimensional Loewner matrices instead of the large-scale ($Q\times K$, $Q\geq K$) $n$-dimensional one, therefore drastically taming the curse of dimensionality, i.e. reducing both the computational effort and the memory needs, and, in addition improving accuracy;
• We show that this decomposition achieves variables decoupling; thus connecting the Loewner framework for rational interpolation of multivariate functions and the Kolmogorov Superposition Theorem (KST), restricted to rational functions. The result is the formulation of KST for the special case of rational functions;
• Connections with KAN neural nets follows (detailed in future work).

Main reference

@article{AGPV:2025,
	Author	= {A.C. Antoulas and I-V. Gosea and C. Poussot-Vassal},
	Title 	= {On the {Loewner} framework, the {Kolmogorov} superposition theorem, and the curse of dimensionality},
	Doi 	= {10.1137/24M1656657},
	Journal = {SIAM Review},
	Volume 	= {64},
	Number 	= {4},
	Pages 	= {737-770},
	Year 	= {2025},
	URL     = {https://doi.org/10.1137/24M1656657}, 
}

Additional support material (presentations, reports, ...)

• A.C. Antoulas presentation: BANFF video
• C. Poussot-Vassal presentation: GT Identification video, slides
• C. Poussot-Vassal presentation: regularly updated slides
• Benchmark results and comparison: Tensor-based multivariate function approximation: methods benchmarking and comparison (arXiv:2506.04791)

The "mLF" MATLAB package

The code (+mlf folder) provided in this GitHub page is given for open science purpose. Its principal objective is to accompany the SIAM Review paper / arXiv paper by the authors, thus aims at being educative rather than industry-oriented. Evolutions (numerical improvements) may come with time. Please, cite the reference above if used in your work and do not hesitate to contact us in case of bug of problem when using it. Below we present an example of use, then functions list are given.

Moreover, for more numerically robust and involved implementation and features, we invite reader and users to refer to the MDSPACK library by MOR Digital Systems.

Dependencies

• MATLAB R2023b or later (tested on this version)
• Toolboxes: symbolic_toolbox may be used for some functionalities

A simple MATLAB code example

We provide two examples to test the approach. Refer to SIAM Review paper / arXiv paper for notations and related equations.

First add the path where the +mlf package is.

%addpath("location_of_mlf") % Add the location of the +mlf package

Quick example

The code demo1.m provides a sample where we use mlf.alg1 and mlf.alg2, standing as implementations of Algorithm 1 and 2 in the above referenced paper. We start by chosing a model in the suggested collection, and construct the tensor along interpolation points.

%%% Pick an example from 1 to 40
CAS         = 1
%%% mlf parameters
alg1_tol    = 1e-9; 
alg2_tol    = 1e-9;
%%% Chose model
[H,infoCas] = mlf.examples(CAS)
n           = infoCas.n;
p_c         = infoCas.p_c;
p_r         = infoCas.p_r;
%%% Data tensor/rand
[y,x,dim]   = mlf.make_tab_vec(H,p_c,p_r);
tab         = mlf.vec2mat(y,dim);

Then, one can use Algorithm 1 (direct).

%%% Alg. 1: direct pLoe [A/G/P-V, 2025]
opt = [];
tic;
opt.ord_tol     = alg1_tol; % SVD tolerance
opt.method_null = 'svd0'; % null space method
opt.method      = 'rec'; % full or recursuve method
% opt.ord_obj     = [];
% opt.ord_N       = 10;
% opt.ord_show    = false;
% opt.data_min    = true;
[r_loe1r,i1_r]  = mlf.alg1(tab,p_c,p_r,opt);
toc

Then, one can use Algorithm 2 (iterative).

%%% Alg. 2: iterative pLoe [A/G/P-V, 2025]
opt = [];
tic
opt.tol         = alg2_tol; % iteration tolerance
opt.method_null = 'svd0'; % null space method
opt.max_iter    = 25; % maximal number of iterations
opt.method      = 'rec'; % full or recursuve method
[r_loe2r,i2_r]  = mlf.alg2(tab,p_c,p_r,opt);
toc

Detailed example

Here is a code that describes the detailed steps on how to deploy the cascaded 1-D Loewner null space construction. Code below is demo0.m.

Then compute the barycentic coefficients for a given function H.

%%% Define a multivariate handle function 
n       = 3; % number of variables
%H       = @(s1,s2,s3) (s3/100-1)*(s2-pi/2)*(s1+atan(2*s2)*tanh(5*(s2-pi)))/(s1^2+s3/10*cos(3*s1)+3)/(s2+10);
H       = @(s1,s2,s3) s2^3/(s1^2+2*s1-s3/2+10);
% /!\ The tolerence is an important parameter when the data are generated from an irrational function
tol_ord = 1e-7; 
% Interpolation points (IP) - separate columns and rows (as Section 3, eq. 13-15)
for ii = 1:n
    p_c{ii} = linspace(-6,6,30);
    dx      = abs(p_c{ii}(2)-p_c{ii}(1))/2;
    p_r{ii} = p_c{ii}+dx;
end

%%% Elementary ingredients
% Check that column/row IP are disjoint (Section 3)
ok                  = mlf.check(p_c,p_r);
% Construct tab_n (Section 3)
tab                 = mlf.make_tab(H,p_c,p_r,true);

%%% Estimate order along each variables and select a subset of IP
ord                 = mlf.compute_order(p_c,p_r,tab,tol_ord,[],5,true);
[pc,pr,W,V,tab_red] = mlf.points_selection(p_c,p_r,tab,ord,false);
w                   = mlf.mat2vec(W);

%%% Recursive null space computation (Section 5, Theorem 5.3)
[c,info]            = mlf.loewner_null_rec(pc,pr,tab_red,'svd0');

%%% Curse of dimensionality: flops & memory estimation (Section 5, Theorem 5.5 & Theorem 5.6)
fprintf('FLOPS\n')
fprintf(' * recursive: %d\n',info.nflop)
fprintf(' * full     : %d\n',length(c)^3)
fprintf('MEMORY\n')
fprintf(' * recursive: %d MB\n',max(ord+1)^2/2^20)
fprintf(' * full     : %d MB\n',prod(ord+1)^2/2^20)

Evaluate simplified function and display results

% Along first and second variables 
% Other variables are randomly chosen between bounds
x1      = linspace(min(p_r{1}),max(p_r{1}),40)+rand(1)/10;
x2      = linspace(min(p_r{2}),max(p_r{2}),41)+rand(1)/10;
[X,Y]   = meshgrid(x1,x2);
rnd_p   = [];
if n > 2; rnd_p = mlf.rand(n-2,p_r(3:end),false); end
for ii = 1:numel(x1)
    for jj = 1:numel(x2)
        param                   = [x1(ii) x2(jj) rnd_p];
        paramStr                = regexprep(num2str(param,36),'\s*',',');
        eval(['tab_ref(jj,ii)   = H(' paramStr ');'] )
        tab_app(jj,ii)          = mlf.eval_lagrangian(pc,w,c,param,false);
    end
end
%
figure
subplot(1,2,1); hold on, grid on
surf(X,Y,tab_app,'EdgeColor','none'), hold on
surf(X,Y,tab_ref,'EdgeColor','k','FaceColor','none')
xlabel('$x_1$','Interpreter','latex')
ylabel('$x_2$','Interpreter','latex')
title('Original vs. Approximation','Interpreter','latex')
axis tight, zlim([min(tab_ref(:)) max(tab_ref(:))]), view(-30,40)
subplot(1,2,2); hold on, grid on, axis tight
imagesc(log10(abs(tab_ref-tab_app)/max(abs(tab_ref(:)))),'XData',x1,'YData',x2)
xlabel('$x_1$','Interpreter','latex')
ylabel('$x_2$','Interpreter','latex')
title('{\bf log}(abs. err.)/max.','Interpreter','latex')
colorbar,

Functions description

mlf.check

Function which checks that column p_c ($P_c^{(n)}$) and row p_r ($P_r^{(n)}$) interpolation points sets are disjoints. Note that both p_c and p_r are $n$-dimensional cells, where each p_c{i} (i=1...n) gathers the interpolatoin points along each variables.

Input

• p_c: column interpolation points $P_c^{(n)}$ ($n$-dimensional cell with double)
• p_r: row interpolation points $P_r^{(n)}$ ($n$-dimensional cell with double)

Output

• ok: tag assessing that interpolation points are disjoints (boolean)

Syntax

ok = mlf.check(p_c,p_r)

mlf.make_tab

Function which constructs the $n$-dimensional tensor $\mathbf{tab}_n$, from the column p_c ($P_c^{(n)}$) and row p_r ($P_r^{(n)}$) interpolation point sets. Note that both p_c and p_r are $n$-dimensional cells, where each p_c{i} (i=1...n) gathers the interpolation points along each variables.

Input

• p_c: column interpolation points $P_c^{(n)}$ ($n$-dimensional cell with double)
• p_r: row interpolation points $P_r^{(n)}$ ($n$-dimensional cell with double)
• show: tag used to show the advancement of the tableau construction (boolean)

Output

• tab: $n$-dimensional tensor $\mathbf{tab}_n$ ($n$-dimensional double)

Syntax

tab = mlf.make_tab(H,p_c,p_r,show);

mlf.compute_order

Function which estimates the orders along each variables.

Input

• p_c: column interpolation points $P_c^{(n)}$ ($n$-dimensional cell with double)
• p_r: row interpolation points $P_r^{(n)}$ ($n$-dimensional cell with double)
• tab: $n$-dimensional tensor $\mathbf{tab}_n$ ($n$-dimensional double)
• ord_tol: normalized singular value threshild for order selection (positive real scalar)
• ord_obj: maximal orders tolerated along each variables ($n\times 1$ integer)
• ord_n: maximal single variable Loewner singular value decompositions per variables (integer scalar)
• show: tag used to show the advancement of the tableau construction (boolean)

Output

• ord: $n$-dimensional vector with order along each variables ($n\times 1$ double)

Syntax

ord = mlf.compute_order(p_c,p_r,tab,ord_tol,ord_obj,ord_n,show);

mlf.points_selection

Function which selects a subset of interpolation points, accordingly to the order along each variables.

Input

• p_c: column interpolation points $P_c^{(n)}$ ($n$-dimensional cell with double)
• p_r: row interpolation points $P_r^{(n)}$ ($n$-dimensional cell with double)
• tab: $n$-dimensional tensor $\mathbf{tab}_n$ ($n$-dimensional double)
• ord: $n$-dimensional vector with order along each variables ($n\times 1$ double)
• square: tag used to set wheather row interpolation points have the same dimension as columns or may be greather; true leads to same number while false allows more row than columns (boolean)

Output

• pc: reduced column interpolation points $P_c^{(n)}$ ($n$-dimnensional cell with double)
• pr: reduced row interpolation points $P_r^{(n)}$ ($n$-dimensional cell with double)
• W: tensor corresponding to $P_r^{(n)}$ evaluation ($n$-dimensional cell with double)
• V: tensor corresponding to $P_c^{(n)}$ evaluation ($n$-dimensional cell with double)
• tab_red: reduced $n$-dimensional tensor $\mathbf{tab}_n$ ($n$-dimensional double)

Syntax

[pc,pr,W,V,tab_red] = mlf.points_selection(p_c,p_r,tab,ord,square);

mlf.loewner_null_rec

Function which construct the Loewner matrix null space (barycentric coefficients) using the recursive approach (without constructing the $N\times N$, or $Q\times K$ Loewner matrix).

Input

• pc: reduced column interpolation points $P_c^{(n)}$ ($n$-dimensional cell with double)
• pr: reduced row interpolation points $P_r^{(n)}$ ($n$-dimensional cell with double)
• tab_red: reduced $n$-dimensional tensor $\mathbf{tab}_n$ ($n$-dimensional double)
• method: method used for null space computation (string); for the time being, we suggest to set this variable to svd0

Output

• c: Loewner matrix null space, being the barycentric coefficients ($N \times 1$ double)
• info: information related to the process (e.g. \mathbf{flop} count)

Syntax

[c,info]  = mlf.loewner_null_rec(pc,pr,tab_red,method);

mlf.eval_lagrangian

Function which evaluates the model in Lagrangian form.

Input

• pc: reduced column interpolation points $P_c^{(n)}$ ($n$-dimensional cell with double)
• w: original vectorized data ealuated at $P_c^{(n)}$ ($n$-dimensional cell with double). In practice we have w = mlf.mat2vec(W);.
• c: Loewner matrix null space, being the barycentric coefficients ($N \times 1$ double)
• param: $n$-variable vector where the model has to be evaluated ($n\times 1$ double)
• show: tag used to show the advancement of the model evaluation (boolean)

Output

• val: value of the Lagrangian model (double)

Syntax

val = mlf.eval_lagrangian(pc,w,c,param,false);

mlf.decoupling

Function which computes the decoupling variables, linking the approach to the Kolmogorov Supperposition Theorem. See Theorem 5.4.

Input

• pc: reduced column interpolation points $P_c^{(n)}$ ($n$-dimensional cell with double)
• info: information given as output of mlf.loewner_null_rec

Output

• Var: single variable weights of the KST $\textbf{Bary}({^1 s})$, $\textbf{Bary}({^1 s})$... ($n$-dimensional cell with double)
• Lag: single variable Lagangian basis of the KST $\textbf{Lag}{{^1 s}}$, $\textbf{Lag}{{^2 s}}$ ... ($n$-dimensional cell with double)
• Bary: Barycentric data information, including $\mathbf{c}^{^1 s}$, $\mathbf{c}^{^2 s}$... and $\mathbf{w}^{^1 s}$, $\mathbf{w}^{^1 s}$...

Syntax

[Var,Lag,Bary] = mlf.decoupling(pc,info);

Feedbacks

Please send any comment to C. Poussot-Vassal (charles.poussot-vassal@onera.fr) if you want to report any bug or user experience issues.

Disclaimer

Once again, this deposit consitutes a research code that accompany the paper mentionned above. It is not aimed to be included in any third party software without the consent of the authors. Authors decline responsabilities in case of problem when applying the code.

Notice also that pathological cases may appear. A more advanced code, to deal with practical and theoretical issues/limitations is currently under developpement by the authors.