I have a dataset consisting of pairs $(x_i, y_i)$. I want to determine the function $f$, so
$$ f(x)f(y) = 1 $$
with the constraint that $f(x) \leq x$, $f'(x) \geq 0$ and $f''(x) \geq 0$.
I was thinking that using splines in some way while constraining the parameters should let me find $f$ satisfying the constraints but I'm unsure how to define the minimization problem in the context of these splines. Any ideas?
Data is available here:
"x" "y"
0.8 1.111
0.76 1.163
0.98 0.92
0.66 1.316
0.9 1
0.78 1.136
1.031 0.87
1.042 0.86
0.85 1.053
1.087 0.82
0.83 1.075
1.099 0.81
0.93 0.97
0.4 2
0.34 2.273
1.053 0.85
1.075 0.83
1 0.9
0.89 1.01
0.91 0.99
0.92 0.98
0.95 0.95
0.82 1.087
0.86 1.042
0.88 1.02
0.41 1.961
0.72 1.22
0.96 0.94
0.7 1.25
1.02 0.88
1.111 0.8
0.81 1.099
1.136 0.78
0.94 0.96
1.19 0.74
0.31 2.439
0.39 2.041
1.25 0.7
0.99 0.91
0.87 1.031
0.97 0.93
1.064 0.84
0.44 1.852
0.84 1.064
0.38 2.083
1.163 0.76
0.68 1.282
0.42 1.923
0.33 2.326
0.75 1.176
0.62 1.389
0.77 1.149
0.61 1.408
0.74 1.19
0.51 1.639
0.6 1.429
0.58 1.471
1.176 0.75
1.124 0.79
0.5 1.667
1.01 0.89
0.46 1.786
1.205 0.73
0.65 1.333
0.48 1.724
0.55 1.538
0.54 1.563
0.37 2.128
0.79 1.124
0.45 1.818
1.149 0.77
0.73 1.205
0.3 2.5
1.22 0.72
0.28 2.632
0.71 1.235
0.35 2.222
0.64 1.351
0.53 1.587
0.63 1.37
0.36 2.174
0.49 1.695
0.32 2.381
0.56 1.515
0.59 1.449
0.67 1.299
0.43 1.887
0.25 2.857
0.69 1.266
0.47 1.754
0.52 1.613
Function $y(x)$ looks monotonical and continued. At the same time, there is a noise due to the discretization. The best approach is the least square method.
$$\color{brown}{\textbf{Linear model}}$$
Let $f(x) = ax+b,$ then should be minimized the function $$G(a,b)=\sum_{i=0}^{n-1}\left((ax_i+b)(ay_i+b)-1\right)^2$$ $$=\sum_{i=0}^{n-1}\left(a^2p_i+abs_i)+b^2-1\right)^2,$$ or $$G(a,b)= a^2(W_{20}a^2+2W_{11}ab+W_{02}b^2)+2a(b^2-1)(W_{10}a+W_{01}b)+W_{00}(b^2-1)^2,$$ where $$ p_i=x_iy_i,\quad s_i - x_i +y_i,\\[4pt] \begin{pmatrix} W_{00} \\ W_{01} \\ W_{10} \\ W_{02} \\ W_{11} \\ W_{20} \end{pmatrix}=\sum_{i=0}^{n-1} \begin{pmatrix} 1 \\ s_i \\ p_i \\ s_i^2 \\ p_is_i \\ p_i^2 \end{pmatrix}= \begin{pmatrix} 91\\ 189.336\\ 78.15834\\ 400.527798\\ 161.52697\\ 67.315936\\ \end{pmatrix} $$ Minimizing leads to the solution with the positive $a:$ $$a_0\approx1.069214,\quad b_0\approx 0.006625$$ (see also Wolram Alpha).
The plot with the values of $\mathbf{\color{lightgreen}{x_i}},$ $\mathbf{\color{blue}{y_i}}$ and $\mathbf{\color{red}{f(x_i)f(y_i)}}$ is shown below. The dotted plots are $f(x_i)$ and $f(y_i).$
$$\color{brown}{\textbf{Quadratic model}}$$
Let $f(x) = ax^2+bx+c,$ then should be minimized the function $$G(a,b,c)=\sum_{i=0}^{n-1}\left((ax_i^2+bx_i+c)(ay_i^2+by_i+c)-1\right)^2$$ $$=\sum_{i=0}^{n-1} \left((a^2p_i+abs_i+b^2)p_i+ac(s_i^2-2p_i)+bcs_i+c^2-1)\right)^2$$ $$=\sum_{i=0}^{n-1} \left((a^2p_i+abs_i+b^2-2ac)p_i+((as_i+b)cs_i+c^2-1))\right)^2$$ $$= a^4W_{40}+2a^3bW_{31}+a^2b^2W_{22}+2(b^2-2ac)(a^2W_{30}+abW_{21})+(b^2-2ac)^2W_{20}$$ $$+2a^2(acW_{22}+bcW_{21}+(c^2-1)W_{20}))$$ $$+2ab(acW_{13}+bcW_{12}+(c^2-1)W_{11}))$$ $$+2(b^2-2ac)(acW_{12}+bcW_{11}+(c^2-1)W_{10}))$$ $$+a^2c^2W_{04}+2abc^2W_{03}+b^2c^2W_{02}+2(c^2-1)(acW_{02}+bcW_{01})+(c^2-1)W_{00},$$ where $$ \begin{pmatrix} W_{00} \\ W_{01} \\ W_{02} \\ W_{03} \\ W_{04} \\ W_{10} \\ W_{11} \\ W_{12} \\ W_{13} \\ W_{20} \\ W_{21} \\ W_{22} \\ W_{30} \\ W_{31} \\ W_{40} \\ \end{pmatrix}=\sum_{i=0}^{n-1} \begin{pmatrix} 1 \\ s_i \\ s_i^2 \\ s_i^3 \\ s_i^4 \\ p_i \\ s_ip_i \\ s_i^2p_i \\ s_i^3p_i \\ p_i^2 \\ s_ip_i^2 \\ s_i^2p_i^2 \\ p_i^3 \\ s_ip_i^3 \\ p_i^4 \\ \end{pmatrix}= \begin{pmatrix} 91\\ 189.336\\ 400.527798\\ 863.931994\\ 1905.452187\\ 78.15834\\ 161.52697\\ 338.997886\\ 724.521647\\ 67.315936\\ 138.255465\\ 288.018453\\ 58.12748\\ 118.699206\\ 50.313345\\ \end{pmatrix} $$ Minimizing via MathCAD with the initial point $$a_1=0,\quad b_1=a_1\approx1.069214,\quad c_1=b_0\approx 0.006625$$ leads to the solution $$a_1\approx 0.314097,\quad b_1\approx 0.75514,\quad c_1=0.006678.$$
The plot with the values of $\mathbf{\color{lightgreen}{x_i}},$ $\mathbf{\color{blue}{y_i}}$ and $\mathbf{\color{red}{f(x_i)f(y_i)}}$ is shown below. The dotted plots are $f(x_i)$ and $f(y_i).$ The values of the quadratic residuals are $g_0=0.228178$ for linear model and $g_1 = 0.006671$ for the quadratic one.
$$\color{brown}{\textbf{Exponential models}}$$
Let $f(x) = e^{h(x)},$ then the condition $f(x)f(y)$ transtorms to $$h(x)+h(y)=0.$$ Let $h(x)=Ax^3+Bx^2+Cx+D,$ then can be minimized the function $$G(a,b)=\sum_{i=0}^{n-1}\left(A(x_i^3+y_i^3)+B(x_i^2+y_i^2)+C(x_i+y_i)+2D\right)^2$$ $$=\sum_{i=0}^{n-1}\left(As_i^3+Bs_i^2+Cs_i+2D-(3As_i+2B)p_i\right)^2,$$ or $$G(a,b)= 9a^2W_{22}+12abW_{21}+4b^2W_{20}-2(3a^2W_{14}+5abW_{13}+(2b^2+3ac)W_{12}+(2bc+6ad)W_{11}$$ $$+4bdW_{10}+a^2W_{06}+2abW_{05}+(2ac+b^2)W_{04}+(4ad+2bc)W_{03}+(4bd+c^2)W_{02}+4cdW_{01}+4d^2W_{00},$$ where $$ p_i=x_iy_i,\quad s_i - x_i +y_i,\\[4pt] \begin{pmatrix} W_{00} \\ W_{01} \\ W_{02} \\ W_{03} \\ W_{04} \\ W_{05} \\ W_{06} \\ W_{10} \\ W_{11} \\ W_{12} \\ W_{13} \\ W_{14} \\ W_{20} \\ W_{21} \\ W_{22} \\ \end{pmatrix}=\sum_{i=0}^{n-1} \begin{pmatrix} 1 \\ s_i \\ s_i^2 \\ s_i^3 \\ s_i^4 \\ s_i^5 \\ s_i^6 \\ p_i \\ p_is_i \\ p_is_i^2 \\ p_is_i^3 \\ p_is_i^4 \\ p_i^2 \\ p_i^2s_i \\ p_i^2s_i^2 \\ \end{pmatrix}= \begin{pmatrix} 91\\ 189.336\\ 400.527798\\ 863.931994\\ 1905.452187\\ 4308.022154\\ 10003.30059\\ 78.15834\\ 161.52697\\ 338.997886\\ 724.521647\\ 1581.452913\\ 67.315936\\ 138.255465\\ 288.018453 \end{pmatrix} $$ Minimizing via MathCAD with the initial point $$a_2=b_2=d_2=0,\quad c_2=1$$ leads to the solution $$a_2\approx 0.00886,\quad b_2\approx-0.075331,\quad c_2\approx0.241154,\quad d_2\approx -0.168262.$$
The plot with the values of $\mathbf{\color{lightgreen}{x_i}},$ $\mathbf{\color{blue}{y_i}}$ and $\mathbf{\color{red}{f(x_i)f(y_i)}}$ is shown below. The dotted plots are $f(x_i)$ and $f(y_i).$ The values of the quadratic residuals are $g_1 = 0.006671$ for the quadratic model and $g2=0.000016$ for the exponential one.
$$\color{brown}{\textbf{Conclusions.}}$$
The most effective model is quadratic one. The linear model has low accuracy, and the exponentional one does not use the negative coefficients.