Rational example

Define a test model and interpolation points

Start by defining your model H $$ H(x_1,x_2) = \frac{x_{1}^2+2x_{2}^2}{x_{1}^2+x_{2}^2}. $$ This is done here by using the mlf.examples routine that embeds multiple examples. We also define the interpolation points along each variable. To fit the Loewner framework, the interpolation points are split into column p_c and row p_r sets. Notice here that the above model is rational along $x_1$ and $x_2$.

CAS         = 55;
[H,infoCas] = mlf.examples(CAS);
n           = infoCas.n;
p_c         = infoCas.p_c;
p_r         = infoCas.p_r;

Build tensor

Then, the tensorized data tab are constructed. Here leading to a 2-D tensor, and more specifically $$ \mathcal T_2^{\otimes} \in \mathbb R^{40\times 40} $$

[y,x,dim]   = mlf.make_tab_vec(H,p_c,p_r);
tab         = mlf.vec2mat(y,dim);

Use the mLF (Alg. 1: direct [A/G/P-V, 2025])

Now we are ready to build the approximation. In what follows, the opt.ord_tol set to $10^{-9}$, is used as the threshold for the single variable order detection. Then the opt.ord_show is used to show the single variable normalized singular value drop (it may help to choose the different orders). The opt.method_null option is used to choose the null space computation method (here SVD with normalization with respect to element with larger magnitude). The opt.method is used to select either the recursive ('rec') or full ('full') method. The algorithm computes g, a handle function; and iloe, gathering some informations about the process and notably the interpolation points $(\lambda_1,\lambda_2)$, the evaluation values $w$ and the barycentric weights $c$. More specifically, the rational approximation takes the following form: $$ g(x_{1},x_{2}) = \dfrac{\sum_{j_1=1}^{k_1}\sum_{j_2=1}^{k_2}\dfrac{c(j_1,j_2)w(j_1,j_2)}{(x_{1}-\lambda_{1}(j_1))(x_{2}-\lambda_{2}(j_2))}}{\sum_{j_1=1}^{k_1}\sum_{j_2=1}^{k_2} \dfrac{c(j_1,j_2)}{(x_{1}-\lambda_{1}(j_1))(x_{2}-\lambda_{2}(j_2))}} $$

opt.ord_tol     = 1e-9;   % SVD tolerance
opt.ord_show    = true;   % show order detection
opt.method_null = 'svd0'; % null space method
opt.method      = 'rec';  % full or recursive method
[g,iloe]        = mlf.alg1(tab,p_c,p_r,opt);

Notice that the singular values decay along both variables $x_1$ and $x_2$ is sharp. This indicates that the model generating the data is rational or polynomial and of degree $d_1=2$ and $d_2=2$. This indicates that the dimension of the nulls psace is $N=(2+1)(2+1)=9$

Verify the null space property

Now it is also possible to verify that, $\mathbf c_2$ (given in iloe.c), the null space computed by the recursive approach (i.e. without constructing the 2-D Loewner matrix), is the solution of $$ \mathbb L_2 \mathbf c_2= 0. $$ To this end one can compute the Loewner matrix and verify its kernel as follows.

tabr        = iloe.tabr;
W           = tabr(1:iloe.ord(1)+1,1:iloe.ord(2)+1);
V           = tabr(iloe.ord(1)+2:end,iloe.ord(2)+2:end);
LL          = mlf.loewnerMatrix(iloe.pc,iloe.pr,W,V);
norm(LL*iloe.c)

Remember that using opt.method='rec' we employ the recursive algorithm (faster and less memory greedy), while opt.method='full' computes the large-scale Loewner matrix and solves the large-scale null space problem. For small scale tensors, 'full' is acceptable but when dimensions and size of the tensor increases, 'rec' become a real game changer. Again, recall that both are mathematically identical but, for some special cases (not clearly identified yet), the recursive approach may lead to numerical issues.

Evaluate the resulting approximate

Now we compare the original and approximate functions with 1000 random draw. The result should be close to machine precision.

for i = 1:1e3
    x_try   = rand(1,2);
    err(i)  = abs(H(x_try)-g(x_try));
end
mean(err)

The below figure compares the original H and approximate g functions along $x_1$ and $x_2$. The frames show the value (left) and the mismatch (right).

Got to the KST form

To go furhter, it is also possible to derive the single variable vector functions $\Phi_1(x_1)$ and $\Phi_2(x_2)$, linking the Loewner framework to the KST, leading the following equivalent formulation of the approximant: $$ g(x_1,x_2) = \dfrac{\sum_{\text{rows}}w\cdot \Phi_1(x_1)\cdot \Phi_2(x_2)}{\sum_{\text{rows}}\Phi_1(x_1)\cdot \Phi_2(x_2)} $$

[Bary,Lag,Cx]   = mlf.decoupling(iloe);
PHI1            = Bary{1}.*Lag{1};
PHI2            = Bary{2}.*Lag{2};
num             = simplify(sum(iloe.w.*PHI1.*PHI2));
den             = simplify(sum(PHI1.*PHI2));
vpa(simplify(num/den),3)
vpa(PHI1,3)
vpa(PHI2,3)

Observe that in this case, the function is perfectly recovered. More specifically, the following KST like vector functions are obtained: $$ \Phi_1(x_1) = \left(\begin{array}{c} \frac{1.11}{x_{1}+1.0}\\ \frac{1.11}{x_{1}+1.0}\\ \frac{1.11}{x_{1}+1.0}\\ -\frac{1.06}{x_{1}+0.0526}\\ -\frac{1.06}{x_{1}+0.0526}\\ -\frac{1.06}{x_{1}+0.0526}\\ \frac{1}{x_{1}-1.0}\\ \frac{1}{x_{1}-1.0}\\ \frac{1}{x_{1}-1.0} \end{array}\right) \text{ and } \Phi_2(x_2) = \left(\begin{array}{c} \frac{1.11}{x_{2}+1.0}\\ -\frac{1.06}{x_{2}+0.0526}\\ \frac{1}{x_{2}-1.0}\\ \frac{1.11}{x_{2}+1.0}\\ -\frac{0.0117}{x_{2}+0.0526}\\ \frac{1}{x_{2}-1.0}\\ \frac{1.11}{x_{2}+1.0}\\ -\frac{1.06}{x_{2}+0.0526}\\ \frac{1}{x_{2}-1.0} \end{array}\right) $$