# On Machine Learning Knowledge Representation In The Form Of Partially
Unitary Operator. Knowledge Generalizing Operator

Vladislav Gennadievich Malyshkin<EMAIL_ADDRESS>Ioffe Institute,
Politekhnicheskaya 26, St Petersburg, 194021, Russia

(December, 22, 2022)

###### Abstract

    $Id: KnowledgeRepresentationAsUnitaryOperator.tex,v 1.449 2022/12/22 06:08:07 mal Exp $

A new form of ML knowledge representation with high generalization power is
developed and implemented numerically. Initial _IN_ attributes and _OUT_ class
label are transformed into the corresponding Hilbert spaces by considering
localized wavefunctions. A partially unitary operator optimally converting a
state from _IN_ Hilbert space into _OUT_ Hilbert space is then built from an
optimization problem of transferring maximal possible probability from _IN_ to
_OUT_ , this leads to the formulation of a new algebraic problem. Constructed
Knowledge Generalizing Operator $\mathcal{U}$ can be considered as a _IN_ to
_OUT_ quantum channel; it is a partially unitary rectangular matrix of the
dimension $\dim(\emph{OUT})\times\dim(\emph{IN})$ transforming operators as
$A^{\emph{OUT}}=\mathcal{U}A^{\emph{IN}}\mathcal{U}^{\dagger}$. Whereas only
operator $\mathcal{U}$ projections squared are observable
$\Braket{\emph{OUT}}{\mathcal{U}}{\emph{IN}}^{2}$ (probabilities), the
fundamental equation is formulated for the operator $\mathcal{U}$ itself. This
is the reason of high generalizing power of the approach; the situation is the
same as for the Schrödinger equation: we can only measure $\psi^{2}$, but the
equation is written for $\psi$ itself.

††preprint: V.M.

##  I Introduction

There are four key elements in any ML approach[1]:

  * •

Attributes selection.

  * •

Knowledge representation.

  * •

Quality criteria (norm).

  * •

Search algorithm to find the solution in knowledge representation space.

Knowledge representation is the most important element as it determines
generalization power of a ML system. The progress in knowledge representation
from linear regression coefficients, perceptron weights[2], statistical
learning[3, 4], and logical approaches[5] to support vector machines[6], rules
and decision trees[7], fuzzy logic[8, 9], and deep learning[10] has been the
direction of ML development within the last two decades.

These approaches, however, share one common feature that limits their
applicability. All of them typically construct a norm, loss function, penalty
function, metric, distance function, etc. on class label (attributes to
predict) difference from the target and perform it’s optimization on training
data. Selection of the norm is a complex task, moreover, the concept of
‘‘norm’’ is of statistical type and cannot be applied in every situation. In
our earlier works[11, 12] we introduced a ‘‘norm-free’’ approach where the
norm was replaced by projection operators. The idea takes inspiration in
quantum mechanics where the outcomes of an observable $f$ (obtained as an
operator’s spectrum $\Ket{f|\psi^{[i]}}=\lambda^{[i]}\Ket{\psi^{[i]}}$) and
the probabilities of outcomes are separated; for a given state $\Ket{\psi}$
the probabilities of $\lambda^{[i]}$ outcomes are obtained as projections to
$\Ket{\psi^{[i]}}$ eigenvectors $\Braket{\psi}{\psi^{[i]}}^{2}$. This approach
comes in two ‘‘flavors’’[13]: interpolatory type (where the outcome is
obtained as regular Radon–Nikodym derivative) and classification type (where
the outcome is obtained as prior weight adjusted Radon–Nikodym derivative, a
‘‘Bayesian’’ style).

While these results are of great interest as they overcome one of the most
difficult problem in ML (norm selection) and produce gauge-invariant
solutions, they, as the other approaches to ML, still have a limitation in
generalization power. The problem with this our approach[13] is that it is
still of ‘‘joint distribution generalization’’ type. Effectively it constructs
a joint distribution of (attributes, class label) pairs and then is trying to
generalize from it. Some ML approaches, such as statistical learning, support
vector machines, rules and decision trees, Bayesian learning, etc. do this
‘‘joint distribution generalization’’ explicitly; the others, such as neural
networks, hidden Markov model, almost all logic models, etc. in fact also do a
‘‘joint distribution generalization’’, but do it implicitly.

The problem with ‘‘joint distribution generalization’’ approaches is that they
can only predict the outcomes that already have corresponding (attributes,
class label) observations in training data. For example if we apply such an
approach to periodic planetary motion – we obtain an accurate prediction, but
applying it to a hyperbolic comet would be a failure as the comet only travel
through the Solar system once. However, both (planet and hyperbolic comet) are
governed by the same Newtonian laws and their motion is the phenomena of the
same kind. Newtonian mechanics has a more powerful generalization than the
‘‘joint distribution generalization’’.

This work is the first work where we go beyond the ‘‘joint distribution
generalization’’ in ML knowledge representation.

##  II Input Data and Simple Models

Whereas the developed approach can be applied to input data of various forms,
for the purpose of comparison with well known models we will be considering
only the data of supervised learning form111 The data can be possibly
“producted” to some order $\mathcal{D}$. For example take $n$ initial $x_{k}$
and construct $x_{\mathbf{k}}=x_{0}^{k_{0}}x_{1}^{k_{1}}\dots
x_{n-1}^{k_{n-1}}$ with multi-index $\mathbf{k}=(k_{0},k_{1},\dots,k_{n-1})$
subject to $\mathcal{D}=\sum\limits_{j=0}^{n-1}k_{j}$. From initial $n$
attributes $x_{k}$ we now obtained
$\mathcal{N}(n,\mathcal{D})=C_{n+\mathcal{D}-1}^{\mathcal{D}}$ attributes
$x_{\mathbf{k}}$ producted to the order $\mathcal{D}$, see [13]. :

$\displaystyle(x_{0},x_{1},\dots,x_{k},\dots,x_{n-1})^{(l)}$
$\displaystyle\to(f_{0},f_{1},\dots,f_{j},\dots,f_{m-1})^{(l)}$ weight
$\omega^{(l)}$ (1) $\displaystyle\mathbf{x}^{(l)}$
$\displaystyle\to\mathbf{f}^{(l)}$

where an attributes vector $\mathbf{x}$ of the dimension $n$ is mapped to a
class label vector $\mathbf{f}$ of the dimension $m$ for all $l=1\dots M$
observations. An average $\Braket{\cdot}$ is defined as the sum over all $M$
observations sample:

$\displaystyle\Braket{1}$ $\displaystyle=\sum\limits_{l=1}^{M}\omega^{(l)}$
(2) $\displaystyle\Braket{h(\mathbf{f})g(\mathbf{x})}$
$\displaystyle=\sum\limits_{l=1}^{M}h(\mathbf{f}^{(l)})g(\mathbf{x}^{(l)})\omega^{(l)}$
(3)

Here $h(\mathbf{f})$ and $g(\mathbf{x})$ are some functions on $\mathbf{f}$
and $\mathbf{x}$, for example a polynomial or Christoffel function $K(x)$ from
(10). In this paper we will be considering the models built on ‘‘moments’’ —
some average of a polynomial function on $x_{k}$ and $f_{j}$; an example of
such an average is $\Braket{x_{k}x_{k^{\prime}}f_{j}f_{j^{\prime}}}$. As a
constant has always to be present in $\mathbf{x}$ and $\mathbf{f}$ bases the
tensor $\Braket{x_{k}x_{k^{\prime}}f_{j}f_{j^{\prime}}}$ includes all lower
order averages such as $\Braket{x_{k}x_{k^{\prime}}}$ and
$\Braket{f_{j}f_{j^{\prime}}}$. Introduce Gram matrices
$G_{kk^{\prime}}^{\mathbf{x}}$ and $G_{jj^{\prime}}^{\mathbf{f}}$ for
$\mathbf{x}$– and $\mathbf{f}$– spaces respectively:

$\displaystyle G_{kk^{\prime}}^{\mathbf{x}}$
$\displaystyle=\Braket{x_{k}x_{k^{\prime}}}$ (4) $\displaystyle
G_{jj^{\prime}}^{\mathbf{f}}$ $\displaystyle=\Braket{f_{j}f_{j^{\prime}}}$ (5)

We will assume that Gram matrices are non–degenerated, otherwise a
regularization to be applied to $\mathbf{x}$ and $\mathbf{f}$ bases, see ‘‘
_Appendix A: Regularization Example_ ’’ of [13].

A few familiar examples. Least squares solution of $\mathbf{f}$ on
$\mathbf{x}$ requires Gram matrix $G_{kk^{\prime}}^{\mathbf{x}}$ and
$\Braket{f_{j}x_{k}}$ moments as input to obtain
$f_{j}(\mathbf{x})=\sum_{k=0}^{n-1}\beta_{k}x_{k}$ as linear system solution:

$\displaystyle\Braket{\left[f_{j}-\sum_{k=0}^{n-1}\beta_{k}x_{k}\right]^{2}}\to\min$
(6) $\displaystyle
f_{j}(\mathbf{x})\approx\sum\limits_{k,k^{\prime}=0}^{n-1}x_{k}G_{kk^{\prime}}^{\mathbf{x};\,-1}\Braket{f_{j}x_{k^{\prime}}}$
(7)

Here $G_{kk^{\prime}}^{\mathbf{x};\,-1}$ is Gram matrix (4) inverse. The (7)
is $m$ different predictors each one is applied to it’s own class label
component $f_{j}$, $j=0\dots m-1$. Least squares knowledge representation
model has limited predictive power and low outlier stability but it is very
easy to implement numerically and obtained solution is gauge-invariant
relatively an arbitrary non–degenerated linear transform of $\mathbf{x}$ and
$\mathbf{f}$:

$\displaystyle x^{\prime}_{k}$
$\displaystyle=\sum\limits_{k^{\prime}=0}^{n-1}T_{kk^{\prime}}x_{k^{\prime}}$
(8a) $\displaystyle f^{\prime}_{j}$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{m-1}T_{jj^{\prime}}f_{j^{\prime}}$
(8b)

This often makes the least squares model the first choice to start data
analysis despite all the drawbacks. The model has the properties similar to
‘‘joint distribution generalization’’ on the support of $\Braket{\cdot}$ and
typically diverges for $\mathbf{x}$ outside of the support interval; it has
low generalization power.

Radon–Nikodym model consists in constructing a weight density
$\psi^{2}_{\mathbf{y}}(\mathbf{x})$ localized at $\mathbf{x}=\mathbf{y}$ and
then averaging $\mathbf{f}$ with it:

$\displaystyle\psi_{\mathbf{y}}(\mathbf{x})$
$\displaystyle=\sqrt{K(\mathbf{y})}\sum\limits_{i,k=0}^{n-1}y_{i}G^{\mathbf{x};\,-1}_{ik}x_{k}=\frac{\sum\limits_{i,k=0}^{n-1}y_{i}G^{\mathbf{x};\,-1}_{ik}x_{k}}{\sqrt{\sum\limits_{i,k=0}^{n-1}y_{i}G^{\mathbf{x};\,-1}_{ik}y_{k}}}=\frac{\sum\limits_{i=0}^{n-1}\psi^{[i]}(\mathbf{y})\psi^{[i]}(\mathbf{x})}{\sqrt{\sum\limits_{i=0}^{n-1}\left[\psi^{[i]}(\mathbf{y})\right]^{2}}}$
(9) $\displaystyle K(\mathbf{x})$
$\displaystyle=\frac{1}{\sum\limits_{i,k=0}^{n-1}x_{i}G^{\mathbf{x};\,-1}_{ik}x_{k}}=\frac{1}{\sum\limits_{i=0}^{n-1}\left[\psi^{[i]}(\mathbf{x})\right]^{2}}$
(10) $\displaystyle f_{j}(\mathbf{x})$
$\displaystyle\approx\frac{\Braket{\psi^{2}_{\mathbf{x}}f_{j}}}{\Braket{\psi^{2}_{\mathbf{x}}}}=\frac{\sum\limits_{i,q,s,k=0}^{n-1}x_{i}G^{\mathbf{x};\,-1}_{iq}\Braket{x_{q}x_{s}f_{j}}G^{\mathbf{x};\,-1}_{sk}x_{k}}{\sum\limits_{i,k=0}^{n-1}x_{i}G^{\mathbf{x};\,-1}_{ik}x_{k}}=\frac{\sum\limits_{i,k=0}^{n-1}\psi^{[i]}(\mathbf{x})\Braket{\psi^{[i]}}{f_{j}}{\psi^{[k]}}\psi^{[k]}(\mathbf{x})}{\sum\limits_{i=0}^{n-1}\left[\psi^{[i]}(\mathbf{x})\right]^{2}}$
(11)

In Eq. (11) the Radon–Nikodym approximation is presented in two bases:
original $x_{k}$, for which $\Braket{x_{i}x_{k}}=G^{\mathbf{x}}_{ik}$, and in
some orthogonalized basis $\Ket{\psi^{[i]}}$ such that
$\Braket{\psi^{[i]}}{\psi^{[k]}}=\delta_{ik}$. Whereas in least squares
approximation (7) the $f_{j}(\mathbf{x})$ is a linear combination of basis
function $x_{k}$, in the Radon–Nikodym approximation (11) it is a ratio of two
quadratic forms on basis function $x_{k}$ with the matrices
$\sum_{q,s=0}^{n-1}G^{\mathbf{x};\,-1}_{iq}\Braket{x_{q}x_{s}f_{j}}G^{\mathbf{x};\,-1}_{sk}$
and $G^{\mathbf{x};\,-1}_{ik}$. By construction it is an averaging with
positive weight222 For a given $\psi$ the normalizing condition is
$1=\Braket{\psi^{2}}$, this is required to properly average an observable
$\Braket{f\psi^{2}}$. In applications, however, the number of “covered”
observations is often also required, for example to estimate possible data
overfitting; the total coverage is $\Braket{1}$ (2). To estimate the number of
observations covered by a given $\psi$ one can use the Christoffel function
$K(\mathbf{x})$ (10) to estimate the coverage as:
$\mathrm{Coverage}_{\psi}\approx\Braket{K\psi^{2}}$. With an expansion of
$K(\mathbf{x})$ in spectrum[13]
$\Ket{K|\psi^{[i]}}=\lambda^{[i]}\Ket{\psi^{[i]}}$ one can obtain an expansion
“by coverage”; this removes the major limitation of the principal components
method: it’s dependence on the scale of $\mathbf{x}$ attributes.
$\Braket{\psi^{2}f_{j}}\big{/}\Braket{\psi^{2}}$ thus the bounds of $f_{j}$
are preserved and the approximation (11) tends to a constant when some
$x_{k}\to\infty$ The calculation requires Gram matrix
$G_{kk^{\prime}}^{\mathbf{x}}$ and $\Braket{x_{k}x_{k^{\prime}}f_{j}}$ moments
as input (compare with $G_{kk^{\prime}}^{\mathbf{x}}$ and
$\Braket{x_{k}f_{j}}$ required for least squares
$f(\mathbf{x})\approx\Braket{\psi_{\mathbf{x}}}{f}\psi_{\mathbf{x}}(\mathbf{x})=\sum_{i,k=0}^{n-1}x_{i}G^{\mathbf{x};\,-1}_{ik}\Braket{x_{k}f}$
approximation); the result is gauge–invariant relatively (8). The (11) is the
solution of ‘‘interpolatory’’ type as it does not take into account ‘‘prior
probabilities’’, see [13] for ‘‘classification’’ type solution with prior
probabilities taken into account, a ‘‘Bayesian style’’.

Figure 1:  A simple demonstration of (9) localized states
$\psi^{2}_{\mathbf{y}}(\mathbf{x})$ for the measure
$\Braket{g}=\int_{-1}^{1}g(x)dx$ and the basis $\mathbf{x}$ constructed from
1D variable $x\in[-1:1]$ as $x_{k}=x^{k}$. The results for the states
localized at $y=\\{-0.6,0,0.4\\}$ are presented as olive, blue, and green
lines respectively. Basis dimension $n$ is chosen as $\\{7,25,50\\}$ for
thick, middle, and thin lines respectively.

A simple demonstration of localized states is presented in Fig. 1. For a
simple chart a multi-dimensional vector $\mathbf{x}$ is constructed from 1D
variable $x\in[-1:1]$ as $x_{k}=x^{k}$. The measure $\Braket{\cdot}$ is taken
as $\Braket{g}=\int_{-1}^{1}g(x)dx$. Then $\psi^{2}_{\mathbf{y}}(\mathbf{x})$
can be considered as a function of scalar $x$ and $y$ as $\mathbf{x}$ and
$\mathbf{y}$ vectors are calculated from the powers of $x$ and $y$. In Fig. 1
we present $\psi^{2}_{-0.6}(x)$, $\psi^{2}_{0}(x)$, and $\psi^{2}_{0.4}(x)$.
As expected the $\psi^{2}_{y}(x)$ density is localized near $x=y$; the
localization becomes stronger with $n$ increase. This chart demonstrates the
main concept behind Radon–Nikodym type of interpolation which is a two–step
process: on the first step a localized state
$\psi^{2}_{\mathbf{y}}(\mathbf{x})$ is built and on the second step the value
of an observable $f$ is evaluated at $\mathbf{y}$ by averaging it with the
weight obtained on the first step:
$f(\mathbf{y})\approx\Braket{\psi^{2}_{\mathbf{y}}f}\big{/}\Braket{\psi^{2}_{\mathbf{y}}}$.
A trivial example of a square wave interpolation using least squares and
Radon–Nikodym is presented in Fig. 2. We see that Radon–Nikodym preserves the
bounds of $f$ and has near interval edge oscillations very much suppressed
because an interpolation of $f$ at $y$ is obtained by averaging $f$ with
always positive weight $d\mu=\psi^{2}_{y}(x)dx$.

Figure 2:  A demonstration of a square wave interpolation (red) by least
squares (blue, Eq. (7)) and Radon–Nikodym (green, Eq. (11)) with the measure
$\Braket{g}=\int_{-1}^{1}g(x)dx$ with $x\in[-1:1]$ for $n=\\{7,25,50\\}$ in
the pictures: left, middle, and right respectively.

###  II.1 Pure Joint Distribution Model

In the section above we considered a simple problem of recovering $\mathbf{f}$
from $\mathbf{x}$ given sampled data (1). The least squares and Radon–Nikodym
estimators (7) and (11) were obtained. They are using individual components of
vector $\mathbf{f}$ as separate class labels; vector class label makes the
study much more difficult than a scalar one. For further development we need,
for attributes $\mathbf{x}$ and class label $\mathbf{f}$ of vector type, to
have estimators of joint distribution $P(\mathbf{x},\mathbf{f})$ probability
and corresponding to it coverage.

There are several possible approaches to unify $\mathbf{x}$ and $\mathbf{f}$.
In [14] the authors introduced a new vector $\mathbf{z}$ of the dimension
$n+m$

$\displaystyle\mathbf{z}$
$\displaystyle=(x_{0},x_{1},\dots,x_{k},\dots,x_{n-1},f_{0},f_{1},\dots,f_{j},\dots,f_{m-1})$
(12)

and constructed Christoffel function from it (this requires all
$\Braket{x_{k}x_{k^{\prime}}}$, $\Braket{f_{j}f_{j^{\prime}}}$, and
$\Braket{x_{k}f_{j}}$ moments). Maximizing Christoffel function on
$\mathbf{f}$ given $\mathbf{x}$ exhibits very promising results. However, a
difficulty with cross-terms arise[13] both in data initial regularization and
in interpretation of the final result.

To deal with vector class label $\mathbf{f}$ and, for further generalization
of Section III below, we will use $\mathbf{f}$-localized states. For sampled
$\mathbf{f}$ data, possibly producted to some order, construct Gram matrix in
$\mathbf{f}$-space (5) and, the same as in (9), build a localized state
$\psi_{\mathbf{g}}(\mathbf{f})$:

$\displaystyle\psi_{\mathbf{g}}(\mathbf{f})$
$\displaystyle=\frac{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}f_{j^{\prime}}}{\sqrt{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}}}=\frac{\sum\limits_{j=0}^{m-1}\psi^{[i]}(\mathbf{g})\psi^{[j]}(\mathbf{f})}{\sqrt{\sum\limits_{j=0}^{m-1}\left[\psi^{[j]}(\mathbf{g})\right]^{2}}}$
(13)

For each observation $l=1\dots M$ consider (14) projection of
$\mathbf{x}^{(l)}$–localized state (9) to $\mathbf{f}^{(l)}$–localized state
(13) then sum it over the entire $l=1\dots M$ sample to obtain the number of
covered observations ${\cal F}^{JDG}$ (note: there is a ‘‘projective’’ factor
$\Braket{f_{j}x_{k}}$ in the expression)

$\displaystyle\Braket{\psi_{\mathbf{g}}}{\psi_{\mathbf{y}}}^{2}$
$\displaystyle=\frac{\left|\sum\limits_{k,k^{\prime}=0}^{n-1}\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j^{\prime}}G^{\mathbf{f};\,-1}_{j^{\prime}j}\Braket{f_{j}x_{k}}G^{\mathbf{x};\,-1}_{kk^{\prime}}y_{k^{\prime}}\right|^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}\sum\limits_{k,k^{\prime}=0}^{n-1}y_{k}G^{\mathbf{x};\,-1}_{kk^{\prime}}y_{k^{\prime}}}$
(14) $\displaystyle\Braket{1}\geq{\cal F}^{JDG}$
$\displaystyle=\sum\limits_{l=1}^{M}\Braket{\psi_{\mathbf{f}^{(l)}}}{\psi_{\mathbf{x}^{(l)}}}^{2}\omega^{(l)}$
(15)

If $\mathbf{x}$ and $\mathbf{f}$ form the same vector space then ${\cal
F}^{JDG}=\Braket{1}$. Otherwise, for example when $\mathbf{x}$ contains the
entire $\mathbf{f}$ plus one more completely random attribute, ${\cal
F}^{JDG}<\Braket{1}$. Since (14) has normalizing terms containing
$G^{\mathbf{f};\,-1}_{jj^{\prime}}$ and $G^{\mathbf{x};\,-1}_{kk^{\prime}}$
matrices in the denominator, to calculate (15) the secondary sampling
technique[15] is required. The Gram matrices $G^{\mathbf{x}}_{kk^{\prime}}$
(4) and $G^{\mathbf{f}}_{jj^{\prime}}$ (5) are calculated first then the
projection (14) is calculated for every observation $l=1\dots M$ and used in
(15) as it were plain observed at observation $l$. Technically this means we
need to calculate the moments of $\mathbf{x}$\- and $\mathbf{f}$\- Christoffel
functions product:
$\Braket{x_{k}f_{j}}{K^{(\mathbf{x})}K^{(\mathbf{f})}}{x_{k^{\prime}}f_{j^{\prime}}}$
(33).

The $\Braket{\psi_{\mathbf{g}}}{\psi_{\mathbf{y}}}^{2}$ can be viewed as joint
distribution of $\mathbf{f}$ and $\mathbf{x}$. For a given $\mathbf{x}$ the
probabilities of various $\mathbf{f}$ can be estimated as

$\displaystyle P(\mathbf{f},\mathbf{x})$
$\displaystyle\approx\Braket{\psi_{\mathbf{f}}}{\psi_{\mathbf{x}}}^{2}$ (16)

The (16) estimates the probability of possible outcome $\mathbf{f}$ given some
fixed value of $\mathbf{x}$; the estimation is based on (attributes, class
label) pairs observed in the training sample. A typical step from here is to
find a subspace of $\mathbf{x}$ providing the best prediction of $\mathbf{f}$,
optimal clustering [13] is a typical approach in this direction. However, we
want to go beyond ‘‘joint distribution generalization’’, beyond finding a
subspace of $\mathbf{x}$ providing the best prediction of $\mathbf{f}$ in
terms of the probability $P(\mathbf{f},\mathbf{x})$ estimated on training
sample. We need a more powerful generalization method, possibly applicable to
not yet observed values of $\mathbf{x}$ and $\mathbf{f}$.

##  III On Knowledge Generalizing Operator

In Section II.1 above we transformed original
$\mathbf{x}^{(l)}\to\mathbf{f}^{(l)}$ data sample (1) to a sequence of
$\mathbf{x}^{(l)}$\- and $\mathbf{f}^{(l)}$\- localized states according to
(9) and (13):

$\displaystyle\psi_{\mathbf{x}^{(l)}}$
$\displaystyle\to\psi_{\mathbf{f}^{(l)}}$ weight $\omega^{(l)}$ (17)

As $\psi_{\mathbf{x}}$ / $\psi_{\mathbf{f}}$ are defined by $n$ / $m$
coefficients before $x_{k}$ / $f_{j}$ the (17) is nothing more than a
transform of the original data (1). This is not a regular linear transform of
$\mathbf{x}$ / $\mathbf{f}$ bases, this is a linear transform with
$G^{\mathbf{x};\,-1}_{kk^{\prime}}$ / $G^{\mathbf{f};\,-1}_{jj^{\prime}}$
matrices followed by normalization to $1$ with Christoffel function as in (9)
/ (13).

The purpose of this transform is to obtain the states we can project to each
other or to some other basis; for example:
$\Ket{\psi_{\mathbf{x}^{(l)}}}=\sum_{k=0}^{n-1}\Ket{\psi^{[k]}}\Braket{\psi^{[k]}}{\psi_{\mathbf{x}^{(l)}}}$
with $\Ket{\psi^{[k]}}$ being an orthogonal full basis in $\mathbf{x}$-space,
$1=\Braket{\psi_{\mathbf{x}^{(l)}}}{\psi_{\mathbf{x}^{(l)}}}^{2}$,
$1=\Braket{\psi_{\mathbf{f}^{(l)}}}{\psi_{\mathbf{f}^{(l)}}}^{2}$,
$1\geq\Braket{\psi_{\mathbf{x}^{(l)}}}{\psi_{\mathbf{f}^{(l)}}}^{2}$, etc. The
main result of Section II.1 was to obtain joint $(\mathbf{x},\mathbf{f})$
distribution (16) and then trying to generalize from it.

Consider a different form of generalization. Let
$\Ket{\psi_{\mathbf{x}^{(l)}}}$, before being used in calculation of joint
distribution, is transformed by a unitary operator $\|\mathcal{U}\|$:

$\displaystyle\mathcal{F}$
$\displaystyle=\sum\limits_{l=1}^{M}\Braket{\psi_{\mathbf{f}^{(l)}}}{\mathcal{U}}{\psi_{\mathbf{x}^{(l)}}}^{2}\omega^{(l)}$
(18)

Contrary to (16) this expression is transforming
$\psi_{\mathbf{x}^{(l)}}(\mathbf{x})$ to some other function in
$\mathbf{x}$-space
$\Ket{\psi(\mathbf{x})}=\Ket{\mathcal{U}|\psi_{\mathbf{x}^{(l)}}(\mathbf{x})}$
and only then projecting the result to actual realization
$\psi_{\mathbf{f}^{(l)}}(\mathbf{f})$ in $\mathbf{f}$-space. In some sense the
$\|\mathcal{U}\|$ can be viewed as a Scattering Amplitude Matrix, as it
relates the _IN_ state $\Ket{\psi_{\mathbf{x}^{(l)}}}$ with the _OUT_ state
$\Ket{\psi_{\mathbf{f}^{(l)}}}$. All the information about what combinations
of attributes $x_{k}$ to be used for prediction now contains in operator
$\|\mathcal{U}\|$. It is called Knowledge Generalizing Operator. The operator
is unitary (to preserve normalizing)
$1=\Braket{\psi}{\mathcal{U}^{\dagger}}{\mathcal{U}|\psi}$.

$\displaystyle\mathcal{U}^{\dagger}\mathcal{U}$ $\displaystyle=\mathds{1}$
(19)

In our model the knowledge is represented in the form of a unitary operator.
This is a very common form in physics: the dynamics in classical mechanics,
electrodynamics, quantum mechanics can be represented as a sequence of
infinitesimal unitary transformations determined by Hamiltonian (or
Lagrangian) of the system. The (18) is an inverse problem: given (1) data find
unitary operator $\|\mathcal{U}\|$ maximizing (18) coverage subject to (19)
constraints. Whereas the (18) is of fourth order in $\Ket{\psi}$, it is of
second order in $\|\mathcal{U}\|$. The constraints (19) is also of second
order in $\|\mathcal{U}\|$. Thus the problem of finding the Knowledge
Generalizing Operator is a variant of quadratically constrained quadratic
program[16] (QCQP).

Before we go further, let us consider a simplified version of (18) to obtain
$\mathcal{F}$ upper limit for ‘‘projective’’ interpretation of operator
$\|\mathcal{U}\|$. Consider the problem of finding (in $\mathbf{x}$-space)
orthogonal basis $\phi^{[i]}$, a subset of full basis $D\leq n$, that
maximizes $\mathcal{F}^{TOT}$:

$\displaystyle\phi^{[i]}$
$\displaystyle=\sum_{k=0}^{n-1}\alpha^{\phi;[i]}_{k}x_{k}$ $\displaystyle
i=0\dots D-1;\enspace D\leq n$ (20) $\displaystyle\delta_{ii^{\prime}}$
$\displaystyle=\Braket{\phi^{[i]}}{\phi^{[i^{\prime}]}}=\sum_{k,k^{\prime}=0}^{n-1}\alpha^{\phi;[i]}_{k}G^{\mathbf{x}}_{kk^{\prime}}\alpha^{\phi;[i^{\prime}]}_{k^{\prime}}$
(21) $\displaystyle\mathcal{F}^{TOT}$
$\displaystyle=\sum\limits_{l=1}^{M}\sum\limits_{i=0}^{D-1}\Braket{\psi_{\mathbf{f}^{(l)}}}{\phi^{[i]}}^{2}\omega^{(l)}$
(22)

Substituting (20) to (13) obtain:

$\displaystyle K^{(\mathbf{f})}(\mathbf{g})$
$\displaystyle=\frac{1}{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}}$
(23) $\displaystyle\Braket{f_{t}}{K^{(\mathbf{f})}}{f_{s}}$
$\displaystyle=\sum\limits_{l=1}^{M}\frac{f^{(l)}_{t}f^{(l)}_{s}}{\sum\limits_{j,j^{\prime}=0}^{m-1}f^{(l)}_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}f^{(l)}_{j^{\prime}}}\omega^{(l)}$
(24) $\displaystyle K^{(\mathbf{f\to x})}_{ik}$
$\displaystyle=\sum\limits_{k^{\prime},t^{\prime},s^{\prime},j^{\prime}=0}^{m-1}\Braket{x_{i}f_{k^{\prime}}}G^{\mathbf{f};\,-1}_{k^{\prime}t^{\prime}}\Braket{f_{t^{\prime}}}{K^{(\mathbf{f})}}{f_{s^{\prime}}}G^{\mathbf{f};\,-1}_{s^{\prime}j^{\prime}}\Braket{x_{k}f_{j^{\prime}}}$
$\displaystyle i,k=0\dots n-1$ (25)

The (24) is $\mathbf{f}$–Christoffel function (23) moments333 One can also
consider $\Braket{\frac{\partial R}{\partial f_{t}}\frac{\partial R}{\partial
f_{s}}}$ with
$R(\mathbf{f})=1/\sqrt{K(\mathbf{f})}=\sqrt{\sum\limits_{j,j^{\prime}=0}^{m-1}f_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}f_{j^{\prime}}}$.
. The (25) allows to present (22) in the form:

$\displaystyle\mathcal{F}^{TOT}$
$\displaystyle=\sum\limits_{i=0}^{D-1}\sum_{k,k^{\prime}=0}^{n-1}\alpha^{\phi;[i]}_{k}K^{(\mathbf{f\to
x})}_{kk^{\prime}}\alpha^{\phi;[i]}_{k^{\prime}}$ (26)

From which we can spectrally expand the $\mathcal{F}^{TOT}$ by solving a
generalized eigenvalue problem with the matrices $K^{(\mathbf{f\to
x})}_{kk^{\prime}}$ and $G^{\mathbf{x}}_{kk^{\prime}}$ in left- and right-
hand sides:

$\displaystyle\sum\limits_{k^{\prime}=0}^{n-1}K^{(\mathbf{f\to
x})}_{kk^{\prime}}\alpha^{\phi;[i]}_{k^{\prime}}$
$\displaystyle=\lambda^{[i]}\sum\limits_{k^{\prime}=0}^{n-1}G^{\mathbf{x}}_{kk^{\prime}}\alpha^{\phi;[i]}_{k^{\prime}}$
(27) $\displaystyle\mathcal{F}^{TOT}$
$\displaystyle=\sum\limits_{i=0}^{D-1}\lambda^{[i]}$ (28)

The (28) is a spectral decomposition of (22), it has at most $m$ non–zero
eigenvalues (the rank of (25) is $m$ or lower, we also assume $m\leq n$). If
$\mathbf{f}$ belongs to a subspace of $\mathbf{x}$ then the sum of these $m$
eigenvalues in (28) is equal to $\Braket{1}$. The (28) takes all possible
vectors from $\mathbf{x}$-space and project them to all
$\Ket{\psi_{\mathbf{f}^{(l)}}}$ summing the coverage, this operation does not
make any inference, it estimates the coverage (18) upper limit for any
norm–preserving projective transform[17], such as $\Ket{\mathcal{U}|\psi}$
(53) or, more generally, (54). The estimation can be obtained from
$K^{(\mathbf{f\to x})}_{kk^{\prime}}$ and $G^{\mathbf{x}}_{kk^{\prime}}$
matrices even without solving the eigenvalue problem (27):

$\displaystyle\mathcal{F}^{TOT}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}K^{(\mathbf{f\to
x})}_{kk^{\prime}}G^{\mathbf{x};\,-1}_{k^{\prime}k}$ (29)

For calculation see com/polytechnik/utils/KGOSolutionVectorXVectorF.java:FTOT
which is used in unit tests.

A simpler approach to construct contributing to coverage subspace
$\Ket{\phi^{[i]}}$ is to notice that in (18) there are scalar products
$\Braket{f_{j}x_{k}}$ of the vectors from $\mathbf{x}$ and $\mathbf{f}$
spaces. Thus we can project the $\mathbf{f}$-space to $\mathbf{x}$-space; to
split $\mathbf{x}$ into two subspaces: $\Ket{\phi^{P;[j]}}$ ‘‘projected’’ (of
the dimension $D\leq m$) and $\Ket{\phi^{O;[k]}}$ ‘‘orthogonal’’ to
$\mathbf{f}$ (of the dimension $n-D$), all vectors from the second one have
zero scalar product with a state in $\mathbf{f}$-space
$\Braket{\phi^{O;[k]}}{\psi_{\mathbf{f}}}=0$, thus the $\Ket{\phi^{O;[k]}}$
does not contribute to coverage (18). For this reason it is sufficient to
consider operator $\|\mathcal{U}\|$ to have the dimension $D\times n$
converting a vector from $\mathbf{x}$-space to $\Ket{\phi^{P;[j]}}$, i.e. to
use $\Ket{\phi^{[i]}}=\Ket{\phi^{P;[i]}}$ as contributing subspace, see
com/polytechnik/utils/TestKGO.java:orthogonalizeU for an implementation.

Let us define operator $\|\mathcal{U}\|$ to be a matrix (in this paper
$u_{sk}$ is considered to be a real matrix, a generalization to a complex
matrix is straightforward) of $D\times n$, $D\leq m$, such that:

$\displaystyle\Ket{\mathcal{U}|x_{k}}$
$\displaystyle=\sum\limits_{s=0}^{D-1}\Ket{\phi^{[s]}}u_{sk}$ (30)

Then (note: there is a ‘‘projective’’ factor $\Braket{f_{j}x_{k}}$ in the
expression, from $\Braket{f_{j}\phi^{[s]}}$)

$\displaystyle\Braket{\psi_{\mathbf{g}}}{\mathcal{U}}{\psi_{\mathbf{y}}}^{2}$
$\displaystyle=\frac{\left|\sum\limits_{k,k^{\prime}=0}^{n-1}\sum\limits_{j,j^{\prime}=0}^{m-1}\sum\limits_{s=0}^{D-1}g_{j^{\prime}}G^{\mathbf{f};\,-1}_{j^{\prime}j}\Braket{f_{j}\phi^{[s]}}u_{sk}G^{\mathbf{x};\,-1}_{kk^{\prime}}y_{k^{\prime}}\right|^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}\sum\limits_{k,k^{\prime}=0}^{n-1}y_{k}G^{\mathbf{x};\,-1}_{kk^{\prime}}y_{k^{\prime}}}$
(31) $\displaystyle\Braket{\phi^{[s]}}{\phi^{[q]}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}u_{sk}\Braket{x_{k}}{x_{k^{\prime}}}u_{qk^{\prime}}$
$\displaystyle s,q=0\dots D-1$ (32)

The (31) is actually (14) but $\Ket{\psi_{\mathbf{y}}}$ is replaced by
$\Ket{\mathcal{U}|\psi_{\mathbf{y}}}$. This is the central concept of
knowledge generalizing operator: the state the inference is based on
$\Ket{\psi_{\mathbf{y}}}$ is transformed by the operator $\|\mathcal{U}\|$
before coupling with the state $\Ket{\psi_{\mathbf{f}}}$ we are looking an
inference to. Partial unitarity constraint (32) corresponds to the fact that
only subspace of the dimension $D\leq m$ can possibly contribute to the
coverage (18). When only a subspace of $\mathbf{x}$ contributes to (35) the
problem to find a unitary matrix $u_{jk}$ becomes highly degenerative. While
the algorithm described in the Appendix A below works well with such a
degenerative problem, it is beneficial for both: computational complexity and
simplicity of result’s interpretation to make the problem less degenerative.
Consider a $\Ket{\phi^{[i]}}$ $i=0\dots D-1$ subspace of the dimension $D\leq
m$. Let us split considered above unitary operator $\|\mathcal{U}\|$ into
$\|\mathcal{U}\|=\|\mathcal{U}^{P}\|+\|\mathcal{U}^{O}\|$ such that
$\|\mathcal{U}^{P}\|$ transforms any $\mathbf{x}$-vector to $\Ket{\phi^{[i]}}$
subspace, and $\|\mathcal{U}^{O}\|$ transforms any $\mathbf{x}$-vector to a
subspace orthogonal to $\Ket{\phi^{[i]}}$ (this split is most easy to perform
if to convert original $\Ket{x_{k}}$ space into direct sum of
$\Ket{\phi^{[i]}}$ and orthogonal to $\Ket{\phi^{[i]}}$ subspaces). Then,
because $0=\Braket{\psi_{\mathbf{f}}}{\mathcal{U}^{O}}{\psi_{\mathbf{x}}}$ for
any $\mathbf{f}$ and $\mathbf{x}$, optimization result of $\mathcal{F}$ does
not depend on $\|\mathcal{U}^{O}\|$, thus it is sufficient to find an operator
$\|\mathcal{U}^{P}\|$ of the dimension $D\times n$ subject to (32) constraint.

To calculate (18) it is convenient to introduce the moments of Christoffel
functions product:

$\displaystyle\Braket{x_{k}f_{j}}{K^{(\mathbf{x})}K^{(\mathbf{f})}}{x_{k^{\prime}}f_{j^{\prime}}}$
$\displaystyle=\sum\limits_{l=0}^{M}\omega^{(l)}\frac{x^{(l)}_{k}x^{(l)}_{k^{\prime}}}{\sum\limits_{q,q^{\prime}=0}^{n-1}x^{(l)}_{q}G^{\mathbf{x};\,-1}_{qq^{\prime}}x^{(l)}_{q^{\prime}}}\cdot\frac{f^{(l)}_{j}f^{(l)}_{j^{\prime}}}{\sum\limits_{s,s^{\prime}=0}^{m-1}f^{(l)}_{s}G^{\mathbf{f};\,-1}_{ss^{\prime}}f^{(l)}_{s^{\prime}}}$
(33)

to write $\mathcal{F}$ in the form (note: there is a ‘‘projective’’ factor
$\Braket{f_{j}x_{k}}$ in the expression, from $\Braket{f_{j}\phi^{[s]}}$)

$\displaystyle S_{sk;s^{\prime}k^{\prime}}$
$\displaystyle=\sum\limits_{j,j^{\prime},q,q^{\prime}=0}^{m-1}\sum\limits_{t,t^{\prime}=0}^{n-1}\Braket{x_{t}f_{j^{\prime}}}{K^{(\mathbf{x})}K^{(\mathbf{f})}}{x_{t^{\prime}}f_{q^{\prime}}}G^{\mathbf{f};\,-1}_{j^{\prime}j}\Braket{f_{j}\phi^{[s]}}G^{\mathbf{x};\,-1}_{kt}G^{\mathbf{f};\,-1}_{q^{\prime}q}\Braket{f_{q}\phi^{[s^{\prime}]}}G^{\mathbf{x};\,-1}_{k^{\prime}t^{\prime}}$
(34) $\displaystyle\mathcal{F}$
$\displaystyle=\sum\limits_{s,s^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{sk}S_{sk;s^{\prime}k^{\prime}}u_{s^{\prime}k^{\prime}}=\sum\limits_{l=1}^{M}\Braket{\psi_{\mathbf{f}^{(l)}}}{\mathcal{U}}{\psi_{\mathbf{x}^{(l)}}}^{2}\omega^{(l)}\xrightarrow[\mathcal{U}]{\quad}\max$
(35)

The $\mathcal{F}$ is a quadratic function on $u_{sk}$; the expression for
$S_{sk;s^{\prime}k^{\prime}}$ can be greatly simplified if $\mathbf{x}$\- and
$\mathbf{f}$\- bases are initially regularized (see [13], ‘‘ _Appendix A:
Regularization Example_ ’’). In an orthogonal basis Gram matrix is a unit
matrix, thus the $G^{\mathbf{x};\,-1}$ and $G^{\mathbf{f};\,-1}$ get removed
in (34).

###  III.1 On Knowledge Generalizing Operator With Different _IN_ and _OUT_
Spaces

In the section above we considered operator $\|\mathcal{U}\|$ as
$\mathbf{x}\to\mathbf{x}$ transform. In (18) the
$\Braket{\psi_{\mathbf{f}}}{\mathcal{U}}{\psi_{\mathbf{x}}}$ was understood as
$\mathbf{x}\to\mathbf{x}$ transform $\Ket{\mathcal{U}|\psi_{\mathbf{x}}}$
followed by projection of the result to $\Ket{\psi_{\mathbf{f}}}$-space (31);
similar ‘‘projective’’ interpretation was used in (14). This interpretation of
$\|\mathcal{U}\|$ lead us to ‘‘contributing subspace’’ $\Ket{\phi^{[s]}}$ (20)
(which is a subspace of $\mathbf{x}$), equation (31) for
$\Braket{\psi_{\mathbf{f}}}{\mathcal{U}}{\psi_{\mathbf{x}}}^{2}$ (it has
$\Braket{f_{j}x_{k}}$ projective factors) and (32) constraints with the
meaning of scalar product invariance. Optimization problem (35) for $u_{jk}$
matrix of the dimension $D\times n$ allows to determine partially unitary
operator $\|\mathcal{U}\|$. This operator has both _IN_ and _OUT_ subspaces
being a subspace of $\mathbf{x}$.

A natural generalization is to consider an operator $\|\mathcal{U}\|$ with
different subspaces for _IN_ and _OUT_ , this way we can avoid any kind of
‘‘projection’’ what would greatly increase generalizing power of the approach.
Let us consider $\mathbf{x}\to\mathbf{f}$ transform directly. Now $u_{jk}$ is
a $m\times n$ matrix transforming a vector from $\mathbf{x}$-space to
$\mathbf{f}$-space

$\displaystyle f_{j}$ $\displaystyle=\sum\limits_{k=0}^{n-1}u_{jk}x_{k}$
$\displaystyle j=0\dots m-1$ (36)

In a common ‘‘projective’’ paradigm the (36) is multiplied by
$x_{k^{\prime}}$, then after taking the average — least squares (7) are
obtained. Now it is different — we cannot take scalar products
$\Braket{f_{j}x_{k}}$ as $\mathbf{f}$ and $\mathbf{x}$ belong to different
Hilbert spaces. We multiply (36) by itself and take the average — obtain (38)
constraint. Substituting (36) to localized state (13) obtain

$\displaystyle\Braket{\psi_{\mathbf{g}}}{\mathcal{U}}{\psi_{\mathbf{y}}}^{2}$
$\displaystyle=\frac{\left|\sum\limits_{k=0}^{n-1}\sum\limits_{j,s=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{js}u_{sk}y_{k}\right|^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}\sum\limits_{k,k^{\prime}=0}^{n-1}y_{k}G^{\mathbf{x};\,-1}_{kk^{\prime}}y_{k^{\prime}}}$
(37) $\displaystyle\Braket{f_{j}}{f_{j^{\prime}}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}\Braket{x_{k}}{x_{k^{\prime}}}u^{*}_{j^{\prime}k^{\prime}}$
$\displaystyle j,j^{\prime}=0\dots m-1$ (38)

thus the optimization problem does not contain any ‘‘projective’’ factors
$\Braket{f_{j}x_{k}}$

$\displaystyle S_{sk;s^{\prime}k^{\prime}}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{m-1}\Braket{f_{j}x_{k}}{K^{(\mathbf{x})}K^{(\mathbf{f})}}{f_{j^{\prime}}x_{k^{\prime}}}G^{\mathbf{f};\,-1}_{js}G^{\mathbf{f};\,-1}_{j^{\prime}s^{\prime}}$
(39) $\displaystyle\mathcal{F}$
$\displaystyle=\sum\limits_{s,s^{\prime}=0}^{m-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{sk}S_{sk;s^{\prime}k^{\prime}}u^{*}_{s^{\prime}k^{\prime}}=\sum\limits_{l=1}^{M}\Braket{\psi_{\mathbf{f}^{(l)}}}{\mathcal{U}}{\psi_{\mathbf{x}^{(l)}}}^{2}\omega^{(l)}\xrightarrow[\mathcal{U}]{\quad}\max$
(40)

This is the equation. The
$\Braket{\psi_{\mathbf{f}}}{\mathcal{U}}{\psi_{\mathbf{x}}}$ is interpreted as
operator $\|\mathcal{U}\|$ relating the states from two different Hilbert
space, a type of memoryless quantum channel, a map between two spaces of
operators. Every admissible transformation $u_{jk}$ must satisfy Gram matrix
invariance condition (38). This condition can be satisfied only for $m\leq n$
since $\Braket{f_{j}}{f_{j^{\prime}}}$ has the rank $m$ and the matrix in the
right hand side has the rank not greater than $n$; in case $m>n$ one can
consider (41) and obtain (43)

$\displaystyle x_{k}$ $\displaystyle=\sum\limits_{j=0}^{m-1}u_{kj}f_{j}$ (41)
$\displaystyle\Braket{x_{k}}{x_{k^{\prime}}}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{m-1}u_{kj}\Braket{f_{j}}{f_{j^{\prime}}}u_{k^{\prime}j^{\prime}}^{*}$
$\displaystyle k,k^{\prime}=0\dots n-1$ (42)
$\displaystyle\Braket{\psi_{\mathbf{g}}}{\mathcal{U}}{\psi_{\mathbf{y}}}^{2}$
$\displaystyle=\frac{\left|\sum\limits_{j=0}^{m-1}\sum\limits_{k,q=0}^{n-1}y_{k}G^{\mathbf{x};\,-1}_{kq}u_{qj}g_{j}\right|^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}\sum\limits_{k,k^{\prime}=0}^{n-1}y_{k}G^{\mathbf{x};\,-1}_{kk^{\prime}}y_{k^{\prime}}}$
(43)

Thus it is sufficient just to swap $\mathbf{x}$ and $\mathbf{f}$ in numerical
calculations. When working in orthogonal bases
$\delta_{kk^{\prime}}=\Braket{x_{k}}{x_{k^{\prime}}}$ and
$\delta_{jj^{\prime}}=\Braket{f_{j}}{f_{j^{\prime}}}$ the matrix elements of
$S_{sk;s^{\prime}k^{\prime}}$ are (33). Also see Appendix B below for possible
adjustment of probability normalizing.

Mapping an operator $A$ between $\mathbf{x}$\- and $\mathbf{f}$\- spaces is
the same transformation
$A^{\mathbf{f}}_{jj^{\prime}}=\sum_{k,k^{\prime}=0}^{n-1}u_{jk}A^{\mathbf{x}}_{kk^{\prime}}u^{*}_{j^{\prime}k^{\prime}}$
as for Gram matrix (38). The optimization problem (40) has the meaning of
finding a quantum channel conveying the highest possible probability from
$\mathbf{x}$–space to $\mathbf{f}$–space. A remarkable feature of this problem
is that it does not contain any $\Braket{f_{j}x_{k}}$ averages! All the
$\mathbf{x}\to\mathbf{f}$ inference (communication between two ends of quantum
channel) now contains only in operator $\|\mathcal{U}\|$ — a matrix $u_{jk}$
of the dimension $m\times n$ to find from optimization problem (40). This is
an important new result. In [13] coverage optimization problem was always
formulated with some kind of $\mathbf{x}\to\mathbf{f}$ projection; if a model
has $\Braket{f_{j}x_{k}}$ terms – it is of ‘‘projective’’ type such as (14),
(31) or (34) above. The (37) and (122) probabilities do not have
$\Braket{f_{j}x_{k}}$ terms; operator $\|\mathcal{U}\|$ directly (36) relates
$\mathbf{x}$\- and $\mathbf{f}$\- spaces subject to (38) scalar product
invariance; it is the only link between _IN_ and _OUT_ spaces. Familiar least
squares expansion (7) satisfies the required constraints (38)

$\displaystyle\Braket{f_{j}}{f_{j^{\prime}}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}\Braket{f_{j}x_{k}}G^{\mathbf{x};\,-1}_{kk^{\prime}}\Braket{x_{k^{\prime}}f_{j^{\prime}}}$
(44)

only when $\mathbf{f}$ is a subspace of $\mathbf{x}$; Proof: select some
orthogonal bases such as $\delta_{kk^{\prime}}=\Braket{x_{k}}{x_{k^{\prime}}}$
and $\delta_{jj^{\prime}}=\Braket{f_{j}}{f_{j^{\prime}}}$, obtain
$1=\sum_{k=0}^{n-1}\Braket{f_{j}x_{k}}^{2}$, i.e. only when
$\mathbf{x}\to\mathbf{f}$ least squares mapping is exact. Note that one can
always apply Appendix (A.5) method of singular values adjustment to obtain a
partially unitary transform from the least squares or any other mapping that
initially does not satisfy the partial unitarity constraints (38).

The optimization considered above has the objective function quadratic on
partially unitary operator $u_{jk}$. There are other objective functions that
are quadratic on partially unitary operator $u_{jk}$ hence all the
optimization above can be applied to them as well. With (36) definition one
can consider it not as probability amplitude mapping
$\psi_{\mathbf{x}}\to\psi_{\mathbf{f}}$, but as plain value mapping
$\mathbf{x}\to\mathbf{f}$. This is essentially (37) without a denominator.
Consider reproducing kernel
$\sum_{j,j^{\prime}=0}^{m-1}f_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}$,
it has a maximum at $\mathbf{f}=\mathbf{g}$, assume $\mathbf{g}$ is taken from
(36), and sum it squared; obtain

$\displaystyle\mathcal{F}$
$\displaystyle=\sum\limits_{s,s^{\prime}=0}^{m-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{sk}S_{sk;s^{\prime}k^{\prime}}u_{s^{\prime}k^{\prime}}=\sum\limits_{l=1}^{M}\left\lgroup\sum\limits_{j,j^{\prime}=0}^{m-1}f_{j}^{(l)}G^{\mathbf{f};\,-1}_{jj^{\prime}}\sum\limits_{k=0}^{n-1}u_{j^{\prime}k^{\prime}}x_{k^{\prime}}^{(l)}\right\rgroup^{2}\omega^{(l)}\xrightarrow[u]{\quad}\max$
(45) $\displaystyle S_{sk;s^{\prime}k^{\prime}}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{m-1}\Braket{f_{j}x_{k}f_{j^{\prime}}x_{k^{\prime}}}G^{\mathbf{f};\,-1}_{js}G^{\mathbf{f};\,-1}_{j^{\prime}s^{\prime}}$
(46)

This creates a different version of $S_{jk;j^{\prime}k^{\prime}}$, a fourth
order moments–type (46) instead of previously used Christoffel functions
product tensor $S_{sk;s^{\prime}k^{\prime}}$ from (39); an important feature
of (46) is that an application of secondary sampling technique is not required
for it’s calculation.

In this setup the conditions on $\Braket{f_{j}f_{j^{\prime}}}$ and
$\Braket{x_{k}x_{k^{\prime}}}$ are put into the constraints (38) and the
$\Braket{f_{j}x_{k}f_{j^{\prime}}x_{k^{\prime}}}$ is put into the objective
function444 In (45) the scalar product of $f_{j}$ and
$\sum_{k=0}^{n-1}u_{jk}x_{k}$ is squared and then averaged over the sample. In
finding the contributing subspace (20) it is averaged over the sample and then
squared. This means the contributing subspace model assumes the factoring
$\Braket{f_{j}x_{k}f_{j^{\prime}}x_{k^{\prime}}}=\Braket{f_{j}x_{k}}\Braket{f_{j^{\prime}}x_{k^{\prime}}}$.
It is similar to Lebesgue quadratures [18], where interchanging of averaging
and taking square produces new result. . The mapping with this new
$S_{sk;s^{\prime}k^{\prime}}$ maps the values, not the probabilities, but the
values are considered to belong to some vector space. The squared term in (45)
is just a scalar product of two vectors. With (122) normalizing both vectors
be of unit length and the maximal value of the objective function is
$\Braket{1}$. In (45) the vectors do not have this normalizing. One can also
consider a ‘‘partially normalized’’ tensor, the one with only
$K^{(\mathbf{f})}$ term in (33) assuming ‘‘average’’–type normalizing for
$x_{k}$ is due to (38).

$\displaystyle S_{sk;s^{\prime}k^{\prime}}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{m-1}\Braket{f_{j}x_{k}}{K^{(\mathbf{f})}}{f_{j^{\prime}}x_{k^{\prime}}}G^{\mathbf{f};\,-1}_{js}G^{\mathbf{f};\,-1}_{j^{\prime}s^{\prime}}$
(47)

###  III.2 Optimization Problem

The problem of finding the Knowledge Generalizing Operator is now reduced to
maximizing (40) coverage $\mathcal{F}$ (defined by the tensor
$S_{sk;s^{\prime}k^{\prime}}$ of diverse possible forms) subject to (38)
constraints; the meaning of the constraints is to preserve scalar product
(Gram matrix). The result is $u_{jk}$ matrix, $j=0\dots m-1;k=0\dots n-1$.
This operator, given some input state (such as localized state
$\Ket{\psi_{\mathbf{x}}}$), uniquely (within a phase) finds the function in
$\mathbf{f}$-space $\Ket{\mathcal{U}|\psi_{\mathbf{x}}}$ (coefficients
$\alpha_{j}$) that predicts the probability (37) of outcome
$\Ket{\psi_{\mathbf{f}}}$:

$\displaystyle P(\mathbf{f})\Big{|}_{\mathbf{x}}$
$\displaystyle=\Braket{\psi_{\mathbf{f}}}{\mathcal{U}}{\psi_{\mathbf{x}}}^{2}=\frac{\left[\sum\limits_{j=0}^{m-1}\alpha_{j}f_{j}\right]^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}f_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}f_{j^{\prime}}}$
(48)

the $\mathbf{f}$ is equal to the value of the outcome we are interested to
determine the probability of. Given $\mathbf{x}$ the probability of some
outcome $\mathbf{f}$ is a squared linear function on $f_{j}$ multiplied by
Christoffel function.

If, however, not the probability but the value of the outcome is required —
the easiest method to obtain it is to consider all possible $\mathbf{f}$ to
find the maximum555 The probability (48) is invariant with respect to
$f_{j}\to const\cdot f_{j}$ for an arbitrary non–zero $const$. Actual values
of $f_{j}$ are determined using the requirement that the constant has always
to be present in $\mathbf{x}$\- and $\mathbf{f}$\- bases. Since the value of
$f_{j}$ corresponding to this specific index ${j:{f_{j}=const}}$ is always
known (a constant), the actual values of all $f_{j}$ are obtained as $f_{j}/C$
where $C=f_{j:{f_{j}=const}}$; see
com/polytechnik/utils/KGOSolutionVectorXVectorF.java:evaluateAt(double []
xorig).  of (48):

$\displaystyle\mathbf{f}$
$\displaystyle:\max\limits_{\mathbf{f}}P(\mathbf{f})\Big{|}_{\mathbf{x}}=\max\limits_{f_{j}}\frac{\left[\sum\limits_{j=0}^{m-1}\alpha_{j}f_{j}\right]^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}f_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}f_{j^{\prime}}}$
(49)

For 1D class label, where $f_{j}=f^{j}$, the problem is reduced to finding the
roots of a polynomial. In general case the problem can be considered as
generalized eigenvalue problem with the matrices
$\alpha_{j}\alpha_{j^{\prime}}$ (a dyadic product of two vectors) and
$G^{\mathbf{f};\,-1}_{jj^{\prime}}$ in the left- and right- hand sides. It has
a single non-zero eigenvalue (51) (equals to the maximal probability),
corresponding eigenvector (50) gives the most probable outcome $\mathbf{f}$.
The maximal probability of the outcome corresponds to the value

$\displaystyle f^{\max P}_{j}$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{m-1}G^{\mathbf{f}}_{jj^{\prime}}\alpha_{j^{\prime}}$
(50) $\displaystyle P(\mathbf{f}^{\max P})\Big{|}_{\mathbf{x}}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{m-1}\alpha_{j}G^{\mathbf{f}}_{jj^{\prime}}\alpha_{j^{\prime}}$
(51)

The $P(\mathbf{f}^{\max P})\Big{|}_{\mathbf{x}}$ is a certainty of the
outcome, the maximal possible value of (48), a $[0:1]$ bounded function. A
difficulty with this approach is that if $\mathbf{f}$ is constructed from a
scalar function, such as $f_{j}=f^{j}$, this relation may not hold exactly in
the result.

Obtained probability formula (48) is of very general form: a linear function
on $f_{j}$ squared divided by a quadratic form on $f_{j}$. It can be obtained
from many different considerations, the difference between models is in
coefficients $\alpha_{j}$. The simplest solution of this type is a ‘‘direct
projection’’ solution of [13], where we take least squares expansion of
$\Ket{f_{j}}$ in $\Ket{x_{k}}$ (7) and substitute obtained
$\mathbf{f}_{LS}(\mathbf{x})$ as the localization point in (13) to obtain
$\Ket{\psi_{\mathbf{f}_{LS}(\mathbf{x})}}$. This is an example to obtain the
probability of (48) form without quantum channel used.

The problem has remarkable invariance features. Consider (17) mappings
$\psi_{\mathbf{x}^{(l)}}\to\psi_{\mathbf{f}^{(l)}}$, $l=1\dots M$ of
$n$-dimensional vector $\psi_{\mathbf{x}^{(l)}}$ to $m$-dimensional vector
$\Ket{\psi_{\mathbf{f}^{(l)}}}$. The vectors are projected to each other with
operator $\|\mathcal{U}\|$, projection absolute value is then squared and all
summed (40) over the entire sample. The major difference from any observable
value–mapping technique is that if we multiply all $\psi_{\mathbf{x}^{(l)}}$
and $\psi_{\mathbf{f}^{(l)}}$ by random phases $\exp(i\varphi^{(l)})$ the
result will be identical! This is the same as in quantum mechanics: a
wavefunction is defined within a phase, wavefunction absolute value squared
defines the probability, but Schrödinger equation is written for the
wavefunction. Similarly, the knowledge generalizing operator $\|\mathcal{U}\|$
is defined (for complex matrix) within a phase, for real matrix – within a
$\pm 1$ factor, but the probability (48) and coverage (40) are equal to
operator $\|\mathcal{U}\|$ projections squared; individual
$\psi_{\mathbf{x}^{(l)}}$ and $\psi_{\mathbf{f}^{(l)}}$ may have arbitrary
phases.

Optimization problem (40) subject to (38) constraints is a variant of QCQP
problem. It has the form: to find an operator $\|\mathcal{U}\|$ optimally
transforming an _IN_ state $\Ket{\psi_{\mathbf{x}}}$ into an _OUT_ state
$\Ket{\psi_{\mathbf{f}}}$ on (17) data, i.e. the ideology is similar to the
one of S-Matrix. Currently we can solve this optimization problem only
numerically. The problem is similar to an eigenvalue problem, see (71). This
is a new algebraic problem:

$\displaystyle S\mathcal{U}$ $\displaystyle=\lambda\mathcal{U}$ (52)

where $S$ is a Hermitian tensor, ‘‘eigenvector’’ $\|\mathcal{U}\|$ is a
partially unitary $m\times n$ matrix, and ‘‘eigenvalues’’ $\lambda$ is a
Hermitian $m\times m$ matrix; functional (40) extremal value is equal to
$\lambda$ spur (the sum of diagonal elements (73)). The mathematical structure
of this eigenvalue–like problem, an ‘‘eigenoperator’’ problem, requires a
separate study and we hope to obtain important new results soon. Currently —
we have a fast, stable to degeneracy iteration algorithm to find a solution
numerically, see Appendix A below.

Considered model assumes the dynamics is determined by a single unitary
operator, possibly partially unitary. For a $\mathbf{x}$-localized pure state
$\|\rho_{\mathbf{x}}\|=\Ket{\psi_{\mathbf{x}}}\Bra{\psi_{\mathbf{x}}}$ a
unitary operator $\|\mathcal{U}\|$ transforms the density matrix to

$\displaystyle\|\widetilde{\rho}_{\mathbf{x}}\|$
$\displaystyle=\|\mathcal{U}|\rho_{\mathbf{x}}|\mathcal{U}^{\dagger}\|$ (53)

Whereas in quantum mechanics evolution operator $\|\mathcal{U}\|$ corresponds
to the Hamiltonian of the system:
$\mathcal{U}=\exp\left[-i\frac{t}{\hbar}H\right]$, in data analysis knowledge
generalizing operator $\|\mathcal{U}\|$ is obtained from optimization problem
(40) subject to (38) constraint. Quantum evolution of (53) form always
transforms a pure state $\|\rho\|=\Ket{\psi}\Bra{\psi}$ to the pure state
$\|\widetilde{\rho}\|=\Ket{\mathcal{U}|\psi}\Bra{\psi|\mathcal{U}^{\dagger}}$,
and a mixed state $\|\rho\|$ to the mixed state $\|\widetilde{\rho}\|$. In
data analysis there is a common situation when a pure state is transformed
into a mixed state, Markov chain is an example. In this case a more general
form of quantum evolution is required[19]:

$\displaystyle\widetilde{\rho}$ $\displaystyle=\sum\limits_{s}B_{s}\rho
B^{\dagger}_{s}$ (54)

with Kraus operators $B_{s}$ satisfying666 Similarly to (38) Kraus operators
$B_{s}$ can also be considered in a “partially unitary”–style with $b_{s;jk}$
matrices of the dimension $m\times n$ satisfying
$\Braket{f_{j}}{f_{j^{\prime}}}=\sum\limits_{s}\sum\limits_{k,k^{\prime}=0}^{n-1}b_{s;jk}\Braket{x_{k}}{x_{k^{\prime}}}b^{*}_{s;j^{\prime}k^{\prime}}$.
The optimization problem (40) then becomes
$\sum\limits_{s}\sum\limits_{l=1}^{M}\Braket{\psi_{\mathbf{f}^{(l)}}}{B_{s}}{\psi_{\mathbf{x}^{(l)}}}^{2}\omega^{(l)}\xrightarrow[B_{s}]{\quad}\max$.

$\displaystyle\sum\limits_{s}B_{s}^{\dagger}B_{s}$ $\displaystyle=\mathds{1}$
(55)

The data we use in this paper is of pure state to pure state mapping (17). For
other type of input data unitary evolution (53) should be replaced by a more
general form (54); one may think about it as a quantum system evolving with
several Hamiltonians at once $B_{s}=\exp\left[-i\frac{t}{\hbar}H_{s}\right]$,
not as about a system evolving with the Hamiltonian $H=\sum\limits_{s}H_{s}$.
The approach is directly generalizable to e.g. probability distribution to
probability distribution mapping: in this case the observations are not
localized states mapping $\psi_{\mathbf{x}^{(l)}}\to\psi_{\mathbf{f}^{(l)}}$,
but corresponding density matrices mapping
$\|\rho_{\mathbf{x}}^{(l)}\|\to\|\rho_{\mathbf{f}}^{(l)}\|$.

Initial $\mathbf{x}^{(l)}\to\mathbf{f}^{(l)}$ input data (1) was converted to
pure state to pure state mapping
$\psi_{\mathbf{x}^{(l)}}\to\psi_{\mathbf{f}^{(l)}}$ (17) to formulate
optimization problem (40) subject to (38) constraints. It is essential from
methodical point of view to discuss what input moments are required for this
problem (to obtain the tensor $S_{jk;j^{\prime}k^{\prime}}$ (39)) and compare
with other models. This is summarized in the table:

Model |  Tensors Required to Calculate
---|---
Least Squares (7) |  $\Braket{x_{k}x_{k^{\prime}}}$, $\Braket{x_{k}f_{j}}$
Radon-Nikodym (11) |  $\Braket{x_{k}x_{k^{\prime}}}$, $\Braket{x_{k}x_{k^{\prime}}f_{j}}$
$\mathbf{x}$ — $\mathbf{f}$ Christoffel function (12) |  $\Braket{x_{k}x_{k^{\prime}}}$, $\Braket{f_{j}f_{j^{\prime}}}$, $\Braket{x_{k}f_{j}}$
Pure Joint Distribution (14) |  $\Braket{x_{k}x_{k^{\prime}}}$, $\Braket{f_{j}f_{j^{\prime}}}$, $\Braket{x_{k}f_{j}}$, $\Braket{x_{k}f_{j}}{K^{(\mathbf{x})}K^{(\mathbf{f})}}{x_{k^{\prime}}f_{j^{\prime}}}$
Partial Unitarity (KGO) (37) |  $\Braket{x_{k}x_{k^{\prime}}}$, $\Braket{f_{j}f_{j^{\prime}}}$, $\Braket{x_{k}f_{j}}{K^{(\mathbf{x})}K^{(\mathbf{f})}}{x_{k^{\prime}}f_{j^{\prime}}}$
Partial Unitarity (KGO) $K^{(\mathbf{f})}$ (47) |  $\Braket{x_{k}x_{k^{\prime}}}$, $\Braket{f_{j}f_{j^{\prime}}}$, $\Braket{x_{k}f_{j}}{K^{(\mathbf{f})}}{x_{k^{\prime}}f_{j^{\prime}}}$
Partial Unitarity (KGO) adj. (122) |  Beyond moments, no $\Braket{x_{k}f_{j}}$ used.

The major difference — Knowledge Generalizing Operator (KGO) is the only model
that does not require ‘‘projective’’ moments $\Braket{x_{k}f_{j}}$; it
requires Gram matrices (4) and (5) of _IN_ and _OUT_ bases and the moments of
the Christoffel functions product (33). These moments can be obtained with an
application of secondary sampling technique[15]: Gram matrices are built
first; then, for every observation $l=1\dots M$, Christoffel function is
calculated and used as it were plain observed at observation $l$. These
moments777 The (122) KGO model goes “beyond moments”. Even with secondary
sampling it is impossible to build from moments the (18) target functional
with the probability (122). Moreover, this problem is not a QCQP problem.  of
two Christoffel functions product are the input used to formulate the problem
(40). For a Christoffel function in some multi-dimensional vector space
$\bm{r}$ (e.g. $\mathbf{x}$ (10) or $\mathbf{f}$ (23)) with $\Braket{\cdot}$
inner product and non-degenerated Gram matrix
$G_{jj^{\prime}}=\Braket{r_{j}}{r_{j^{\prime}}}$ there is a $1/r^{2}$
asymptotic:

$\displaystyle K(\bm{r})$
$\displaystyle=\frac{1}{\Braket{\bm{r}}{G^{-1}}{\bm{r}}}=\frac{1}{\sum\limits_{j,j^{\prime}}r_{j}G^{-1}_{jj^{\prime}}r_{j^{\prime}}}$
(56) $\displaystyle K(\bm{r})$ $\displaystyle\sim 1/r^{2}$
$\displaystyle\text{for }\enspace r\to\infty$ (57)

The same $1/r^{2}$ long–range interaction presents in Coulomb’s law or
Newton’s law of gravitation. With (57) asymptotic the Christoffel function can
be viewed as a form of ‘‘long–range $1/r^{2}$ interaction’’, an anisotropic
gravity–like law of data analysis. These non-local features, along with
eigenproblem (70) of the dimension $Dn$ and SVD (74) (or Gram matrix
eigenproblem (90)) that are required on every iteration, substantially slow
down the algorithm when implemented without optimization. At this point,
however, the goal is not to build a fast algorithm, but to understand all the
benefits and drawbacks of ML knowledge representation in the form of partially
unitary operator. Let us do a demonstration.

##  IV A Demonstration Of Knowledge Generalizing Operator Application

In this section we are going to present several demonstrations of
$\mathbf{f}(\mathbf{x})$ calculation using (49). The $\mathbf{f}$ and
$\mathbf{x}$ are treated as linear spaces, a basis for wavefunction, with
partially unitary operator $u_{jk}$ mapping (36). The result is invariant
relatively $\mathbf{f}\to C\cdot\mathbf{f}$. To obtain actual value of
$\mathbf{f}$ — it should be normalized to const. The constant has always to be
present in both $\mathbf{f}$– and $\mathbf{x}$– bases. Thus

$\displaystyle\mathbf{f}^{actual}$
$\displaystyle=\frac{\mathbf{f}}{f_{j:{f_{j}=const}}}$ (58)

In this equation the nominator is a linear function on $\mathbf{x}$ (36) and
the denominator, the const–component of $\mathbf{f}$, possibly also is a
linear function on $\mathbf{x}$. Thus the value obtained from partially
unitary operator mapping is a ratio of two linear functions on $\mathbf{x}$.
The least squares (7) always maps a constant to a constant, thus when $u_{jk}$
is a least squares mapping the denominator in (58) is always a constant. In
Radon–Nikodym mapping (11) the nominator is a quadratic form on $\mathbf{x}$
and the denominator is a positive quadratic form on $\mathbf{x}$; the
denominator is never zero. In (58) the nominator and the denominator are both
linear functions on $\mathbf{x}$ of most general form. The divergences coming
from denominator’s zeroes are important new features of the approach. In least
squares – these zeroes are on the infinity. Denominator’s zeroes may come
either from deep internal properties of the model or from sub-optimal solution
of the optimization problem (or badly chosen objective function).

The objective function is determined by the tensor
$S_{jk;j^{\prime}k^{\prime}}$. Whereas properly normalized probability (122)
lead to a non–QCQP problem, the original Christoffel (39), the adjusted number
of degrees of freedom Christoffel (128), $\mathbf{f}$–Christoffel (47), and
plain $\Braket{f_{j}x_{k}f_{j^{\prime}}x_{k^{\prime}}}$ (46) have
$S_{jk;j^{\prime}k^{\prime}}$ tensor readily available and the optimization
problem (65) with the constraints (66) can be formulated and solved
numerically.

Among available $S_{jk;j^{\prime}k^{\prime}}$ versions the
$\mathbf{f}$–Christoffel (47) has the most ‘‘usual’’ properties. For example
the (39) or (128), when run with a data of exact $\mathbf{x}\to\mathbf{f}$
homomorphism can possibly give a higher $\mathcal{F}$ on non–exact mapping due
to unusual localized states normalizing. For this reason all the
demonstrations below will be performed with $\mathbf{f}$–Christoffel
$S_{jk;j^{\prime}k^{\prime}}$ (47).

Figure 3:  For a data with known exact $f=x$ solution, when numerical method
does not find it – it is possible to have zeroes in (58) and corresponding
poles in the behavior.

Consider a trivial mapping with the measure $\Braket{g}=\int_{-1}^{1}g(x)dx$
and the basis $\mathbf{x}$ constructed from 1D variable $x\in[-1:1]$ as
$x_{k}=x^{k}$, $k=0\dots 6;n=7$, and $f_{j}=x_{j}$ for $j=0\dots 4;m=5$. The
solution is trivial: take first $m$ components of $x_{k}$ and regularize; then
use them for both: $\mathbf{x}$ and $\mathbf{f}$. However, when the numerical
algorithm cannot find this exact solution we can observe a deviation from
exact match. In Fig. 3 the exact solution along with two approximate solutions
of different quality are presented. A not very accurate approximate numerical
solution may give poles corresponding to the zeroes in (58) denominator
(clearly observed for $f^{\mathrm{approx}\,2}$ near interval edge).

In Fig. 4 a square wave step function (the same as in Fig. 2) is presented
with the same measure and basis; $n=7$. The $f$ takes only two values since
the only available $m$ is $m=2$. The exact solution was difficult to obtain
numerically as the problem is substantially degenerated. We present three
approximate solutions. The blue line is regular least squares (7). Light blue
is the same least squares mapping (7) adjusted with (75) to partial unitarity.
Green — maximal eigenvalue (67) solution adjusted to partial unitarity with
(75). One can see that partial unitarity adjustment makes little changes to
least squares solution. For adjusted maximal eigenvalue solution the (58)
denominator poles are close to the support of $\mathbf{x}$, this creates two
artifacts in $\mathbf{f}$. Note almost exact $f=1$ matching in the center.

Figure 4: A square wave step function (the same as in Fig. 2), with least
squares (blue), least squares with (75) SVD adjustment (light blue), and
maximal eigenvalue with (75) SVD adjustment (green).

Consider a 2D example. Let us take an image and consider it as a
two–dimensional basis mapping a pixel coordinate $(x,y)$ to gray intensity
$f$.

$\displaystyle(x,y)^{(l)}$ $\displaystyle\to f^{(l)}$ weight $\omega^{(l)}=1$
(59) $\displaystyle\mathbf{k}$ $\displaystyle=(k_{x},k_{y})$ (60)
$\displaystyle x_{\mathbf{k}}$ $\displaystyle=x^{k_{x}}y^{k_{y}}$
$\displaystyle 0\leq k_{x}\leq n_{x}-1;{\quad}0\leq k_{y}\leq n_{y}-1$ (61)
$\displaystyle f_{j}$ $\displaystyle=f^{j}$ $\displaystyle j=0\dots m-1$ (62)

This forms a (1) basis888 For numerical stability it is better to use
argument–scaled Chebyshev polynomials rather than monomials powers
$x^{k_{x}}y^{k_{y}}$ and $f^{j}$. of $n=n_{x}n_{y}$ and $m$ dimensions. Let us
construct an operator $u_{jk}$ mapping $\mathbf{x}\to\mathbf{f}$. A simple
example is least squares (7), it creates a familiar image expansion similar to
Fourier series. However, we are interested in operators $u_{jk}$ satisfying
all partial unitarity constraints (38). A simple variant of
constraint–satisfying operator can be obtained from any $u_{jk}$ operator
applying Appendix A.5 algorithm.

Figure 5:  A demonstration of image interpolation with $n_{x}=n_{y}=25$,
$m=5$. Top row: original image, least squares (7) interpolated, and the same
least squares adjusted with (75) to partial unitarity constraints. Middle row:
optimization (40) with simplified constraints (69) (the state of maximal
eigenvalue), the same one adjusted with (75) to partial unitarity constraints,
and optimization result with Section A.4 algorithm. Bottom row: The
probability (48) is calculated at actual $\mathbf{f}$, white $P=1$, black
$P=0$. It is calculated for: least squares (7) (‘‘direct projection’’ model of
[13]), the state of maximal eigenvalue (unadjusted), and Section A.4
algorithm.

In Fig. 5 (top row) we present original image, least squares expansion and
constraint–adjusted least squares for $n_{x}=n_{y}=25$, $m=5$. The
constraint–adjusted least squares is very similar to the original least
squares. The least squares operator maps pixel coordinates to gray intensity,
not the localized states wavefunction. When an operator is optimized to map
the wavefunctions this may cause poles in values, the zeroes of (58)
denominator. It is trivially to construct a partially unitary operator
$u_{jk}$ preserving the constant: construct a partially unitary operator
mapping $\mathbf{x}$-space without const to $\mathbf{f}$-space without const
$\mathbf{x}\backslash C\to\mathbf{f}\backslash C$, then do a direct sum with
$C\to C$ mapping. We do not perform such a transform specifically to observe
the poles in (58). We present three pictures, corresponding to $u_{jk}$
operators differently optimizing (40) with $\mathbf{f}$–Christoffel tensor
(47). In Fig. 5 (middle row) we present the results corresponding to these
three $u_{jk}$: optimizing (40) with simplified constraints (69), the same one
adjusted with (75) to partial unitarity, and optimization result with Section
A.4 algorithm (overall the best optimization algorithm we have so far). Left
two pictures in the middle row — a simple solutions (based on trivial approach
of maximal eigenvalue state), they have noticeable $1/n_{\\{x,y\\}}$ scale
artifacts. The last one is very close to the global maximum of (40) and
‘‘mixes’’ the modes much stronger . The poles of (58) separate the regions and
the structure of these ‘‘separators’’ can be a subject of our future research.

The developed approach works with probabilities, not with the values. For this
reason it is of interest to present the probability (48) at given known
outcome $\mathbf{f}=\mathbf{f}^{(l)}$. The result is presented in Fig. 5, the
bottom row. The probability is scaled as white $P=1$, black $P=0$. It is
presented in the bottom row for three algorithms: least squares (7), the state
of maximal eigenvalue (unadjusted), and Section A.4 algorithm.

The method to overcome noticeable $1/n_{\\{x,y\\}}$ artifacts in Fig. 5 is to
use properly normalized states (122). In most general form it can be
considered as an unconstrained optimization problem. Given sampled data (1)
find a linear transform $u_{jk}$ (63), a general form matrix of the dimension
$m\times n$, maximizing (64)

$\displaystyle f_{j}$ $\displaystyle=\sum\limits_{k=0}^{n-1}u_{jk}x_{k}$
$\displaystyle j=0\dots m-1$ (63) $\displaystyle\mathcal{F}$
$\displaystyle=\sum\limits_{l=1}^{M}\Braket{\psi_{\mathbf{f}^{(l)}}}{\psi_{u(\mathbf{x}^{(l)})}}^{2}\omega^{(l)}\xrightarrow[u]{\quad}\max$
(64)

Here the $\Ket{\psi_{\mathbf{g}}}$ is the state (13) localized at
$\mathbf{f}=\mathbf{g}$, and $\Ket{\psi_{u(\mathbf{x})}}$ is also
$\mathbf{f}$–localized state (13) with the localization point $\mathbf{g}$
determined by (63) linear mapping. When expanded
$\Braket{\psi_{\mathbf{g}}}{\psi_{u(\mathbf{y})}}^{2}$ is (122). The objective
function (64) is the total probability transferred from $\mathbf{x}$–space to
$\mathbf{f}$–space; this is an unconstrained problem. In this most general
form the problem is not a QCQP problem and it is difficult to solve
numerically; the difficulty is that with $\Ket{\psi_{u(\mathbf{x})}}$ state
the operator $u_{jk}$ enters (through localization point $\mathbf{g}$) both
the nominator and the denominator of (13), what makes the optimization problem
(64) not a QCQP problem. The problem can be substantially simplified when the
$u_{jk}$ mapping is considered to be a partially unitary transform (38) to
obtain a QCQP problem. The problem can be further approximated by splitting
the solution into two steps: selecting the contributing subspace $\phi_{k}$ of
the dimension $m$, then constructing a unitary (not partially unitary) mapping
from the contributing subspace to $f_{j}$. A simple projective approach is
presented above in Eq. (26) or, more generally, in the Appendix C below. A
simple solution of this type is the ‘‘direct projection’’ model of [13] where
the localization point is determined from plain least squares (7) to obtain
the state $\Ket{\psi_{\mathbf{f}_{LS}(\mathbf{x})}}$. The probability
$\Braket{\psi_{\mathbf{f}^{(l)}}}{\psi_{\mathbf{f}_{LS}(\mathbf{x}^{(l)})}}^{2}$
of the ‘‘direct projection’’ model is presented in Fig. 5 (leftmost in the
bottom row).

These demonstrations make us to conclude that partial unitary mapping is a
rich form of knowledge representation with a high generalizing power, however
a more study is required.

##  V Conclusion

The developed knowledge generalizing operator concept is similar to the
S-Matrix approach since it is an operator optimally transforming an _IN_ state
$\Ket{\psi_{\mathbf{x}}}$ into an _OUT_ state $\Ket{\psi_{\mathbf{f}}}$. As
any wavefunction in ML is known within an arbitrary phase the equation for the
operator must include only observable values. The problem we consider is to
recover $\|\mathcal{U}\|$ from all it’s projections squared, from the
probabilities (37). The condition of operator’s optimality is (40) coverage
maximization on (17) data; it is a new kind of algebraic problem (52) — the
equation to determine the $\|\mathcal{U}\|$. The situation is the same as with
the Schrödinger equation: the equation is written for $\psi$, but only
$\psi^{2}$ is observable. This is the difference between our and all other ML
knowledge representation techniques where knowledge representation
characteristics are observable values. If a model relates an initial
observable and the final observable then it is a ‘‘joint distribution model’’,
it cannot predict something that has not been already observed in the training
data. Knowledge generalizing operator relates the amplitude of the initial
state to the amplitude of the final state. This is the very feature that
creates generalization. The same is in quantum mechanics: $\psi^{2}$ vs
$\psi$; whereas a mapping of $\psi^{2}$ is meaningless, the mapping of $\psi$
determines the dynamics of a system.

Considered maximization problem (40) is a simple example of knowledge
generalizing operator technique: for observations $l=1\dots M$ convert
$\mathbf{x}^{(l)}\to\mathbf{f}^{(l)}$ to
$\psi_{\mathbf{x}^{(l)}}\to\psi_{\mathbf{f}^{(l)}}$, then reconstruct
$\|\mathcal{U}\|$ from it’s projections squared
$\Braket{\psi_{\mathbf{f}^{(l)}}}{\mathcal{U}}{\psi_{\mathbf{x}^{(l)}}}^{2}$.
The problem can be generalized by considering, instead of $l$, $\mathbf{x}$,
$\mathbf{f}$, and $\Braket{\cdot}$, the structures generalizing the concepts
of set, vector, and measure. In the most general form it can be formulated as:
for $\psi\in S_{x}$ and $\varphi\in S_{f}$ recover partially unitary operator
$\|\mathcal{U}\|$ from it’s projections squared $\sum\limits_{l\in
M}\omega^{(l)}\Braket{\varphi^{(l)}}{\mathcal{U}}{\psi^{(l)}}^{2}\xrightarrow[\mathcal{U}]{\quad}\max$.
The problem can be further generalized by considering mixed states
$\|\rho\|\in S_{x}$ and $\|\varrho\|\in S_{f}$ and recovering Kraus operators
$B_{s}$ (54) from projections squared: $\sum\limits_{l\in
M}\omega^{(l)}\sum\limits_{s}\mathrm{Spur}\|\varrho^{(l)}|B_{s}|\rho^{(l)}|B_{s}^{\dagger}\|\xrightarrow[B_{s}]{\quad}\max$.

There is another interesting twist to the considered problem of finding a
partially unitary matrix $u_{jk}$ of the dimension
$\dim(\emph{OUT})\times\dim(\emph{IN})$ mapping operators from
$\Ket{\emph{IN}}$ to $\Ket{\emph{OUT}}$. Consider the problem: for
$\dim(\emph{OUT})<\dim(\emph{IN})$ select $\dim(\emph{OUT})$ input attributes
out of all $\dim(\emph{IN})$ available that maximize some correctness
condition which is a function of all selected attributes. For all interesting
correctness conditions this problem is typically a one of NP–complete type.
There is a single correctness function (least squares) that can be trivially
solved. Maximization of total matched probability (18) among all partially
unitary operators $u_{jk}$ also selects $\dim(\emph{OUT})$ inputs from all
$\dim(\emph{IN})$ available. This is a new algebraic problem (52). Found
mapping $u_{jk}$ can be viewed as a solution to attributes selection problem
with correctness conditions somewhere ‘‘in between’’ least squares and
NP–complete, for example there is a simple subspace selection approach (134) —
then a problem of unitary mapping (not partially unitary) can be directly
solved.

##  Appendix A A Numerical Solution to Find the Knowledge Generalizing
Operator

The problem we consider is a QCQP problem to maximize (65) subject to (66)
constraint.

$\displaystyle\mathcal{F}$
$\displaystyle=\sum\limits_{l=1}^{M}\Braket{\psi_{\mathbf{f}^{(l)}}}{\mathcal{U}}{\psi_{\mathbf{x}^{(l)}}}^{2}\omega^{(l)}=\sum\limits_{j,j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}S_{jk;j^{\prime}k^{\prime}}u_{j^{\prime}k^{\prime}}\xrightarrow[u]{\quad}\max$
(65) $\displaystyle\Braket{f_{j}}{f_{j^{\prime}}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}\Braket{x_{k}}{x_{k^{\prime}}}u_{j^{\prime}k^{\prime}}$
$\displaystyle j,j^{\prime}=0\dots D-1$ (66)

Without loss of generality we put
$\delta_{kk^{\prime}}=\Braket{x_{k}}{x_{k^{\prime}}}$ and
$\delta_{jj^{\prime}}=\Braket{f_{j}}{f_{j^{\prime}}}$ as we can always choose
an orthogonal basis by applying, for example, an orthogonalization of
Gram–Schmidt type. Contrary to other methods (e.g. regular principal
components) the result obtained with knowledge generalizing operator is
invariant with respect to (8) transform of input data, thus it does not depend
on initial regularization. The problem becomes:

$\displaystyle\mathcal{F}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}S_{jk;j^{\prime}k^{\prime}}u_{j^{\prime}k^{\prime}}\xrightarrow[u]{\quad}\max$
(67) $\displaystyle\delta_{jj^{\prime}}$
$\displaystyle=\sum\limits_{k=0}^{n-1}u_{jk}u_{j^{\prime}k}$ $\displaystyle
j,j^{\prime}=0\dots D-1$ (68)

Consider the squared Frobenius norm of matrix $u_{jk}$ to be a ‘‘simplified
constraint’’:

$\displaystyle\sum\limits_{j=0}^{D-1}\sum\limits_{k=0}^{n-1}u^{2}_{jk}=D$ (69)

This is a ‘‘partial’’ constraint (it is the sum of all (68) diagonal
elements). For this ‘‘partial’’ constraint optimization problem (67) can be
readily converted to an eigenvalue problem that can be directly solved. The
main idea is to adjust this ‘‘preliminary’’ solution to satisfy the full set
of (68) constraints and then calculate new values of Lagrange multipliers.
Performing several iterations the process possibly converge to (67) maximum
with all the required constraints (68) satisfied. In [13] a similar technique
has been tried for a unitary operator (19). The (68) corresponds to partially
orthogonal operator (partially unitary real matrix): $D\leq n$.

Consider Lagrange multipliers $\lambda_{jj^{\prime}}$, a matrix of $D\times D$
dimension, to approach optimization problem (67) with the constraints (68)

$\displaystyle\sum\limits_{j,j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}S_{jk;j^{\prime}k^{\prime}}u^{*}_{j^{\prime}k^{\prime}}+\sum\limits_{j,j^{\prime}=0}^{D-1}\lambda_{jj^{\prime}}\left[\delta_{jj^{\prime}}-\sum\limits_{k^{\prime}=0}^{n-1}u_{jk^{\prime}}u^{*}_{j^{\prime}k^{\prime}}\right]\xrightarrow[u]{\quad}\max$
(70)

Despite the matrix $u_{jk}$ being real we write it in a ‘‘complex’’ form to
variate separately over $u_{jk}$ and $u^{*}_{jk}$. The tensor
$S_{jk;j^{\prime}k^{\prime}}=S^{*}_{j^{\prime}k^{\prime};jk}$ is Hermitian.
The variations

$\displaystyle 0$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{D-1}\sum\limits_{k^{\prime}=0}^{n-1}u_{j^{\prime}k^{\prime}}S_{j^{\prime}k^{\prime};jk}-\sum\limits_{j^{\prime}=0}^{D-1}\lambda_{j^{\prime}j}u_{j^{\prime}k}$
(71a) $\displaystyle 0$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{D-1}\sum\limits_{k^{\prime}=0}^{n-1}S_{jk;j^{\prime}k^{\prime}}u^{*}_{j^{\prime}k^{\prime}}-\sum\limits_{j^{\prime}=0}^{D-1}\lambda_{jj^{\prime}}u^{*}_{j^{\prime}k}$
(71b)

are consistent only when $\lambda_{jj^{\prime}}$ is a Hermitian matrix

$\displaystyle\lambda_{jj^{\prime}}$ $\displaystyle=\lambda^{*}_{j^{\prime}j}$
(72)

From (71) it follows that the functional (67) extremal value is equal to the
spur of Lagrange multipliers matrix $\lambda_{jj^{\prime}}$:

$\displaystyle\max{\mathcal{F}}$
$\displaystyle=\sum\limits_{j=0}^{D-1}\lambda_{jj}$ (73)

An iteration algorithm finding the maximum of (67) subject to (68) constraints
is:

  1. 1.

Take initial $\lambda_{ij}$ and solve optimization problem (70) with respect
to $u_{jk}$ subject to partial constraint (69). Solution method — an
eigenvalue problem of $Dn$ dimension in a vector space formed by writing all
$u_{jk}$ matrix elements in a vector, row by row. The result: $p=0\dots Dn-1$
eigenvalues ${\mathcal{F}}^{[p]}$ and corresponding matrices $u_{jk}^{[p]}$
reconstructed back from the eigenvectors, row by row.

  2. 2.

To select the $u_{jk}$ among all $Dn$ eigenstates one need to try a number of
them, selecting the ones providing a large value of the original functional.
Taking only the state of the maximal eigenvalue typically gives a local
maximum. Chosen $u_{jk}$ is not partially unitary as the constraint (69) is a
subset of the full ones (68). Expand $u_{jk}$ in SVD:

$\displaystyle u_{jk}$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{D-1}\sum\limits_{k^{\prime}=0}^{n-1}U_{jj^{\prime}}\Sigma_{j^{\prime}k^{\prime}}V^{\dagger}_{k^{\prime}k}$
(74)

and adjust all SVD numbers to $\pm 1$. The $\Sigma_{jk}=\delta_{jk}$ is
typically the best option as this is the minimal change (initial $\Sigma_{jj}$
are positive). Obtained

$\displaystyle\widetilde{u}_{jk}=\sum\limits_{s=0}^{\min(D,n)-1}U_{js}V^{\dagger}_{sk}$
(75)

is a partially unitary matrix satisfying all the constraints (68). This
$\widetilde{u}_{jk}$ becomes the next iteration $u_{jk}$ of the solution.
Because of ${u}_{jk}\to\widetilde{u}_{jk}$ adjustment the value of
${\mathcal{F}}$ becomes less optimal. There are other methods to adjust the
$u_{jk}$ to satisfy the full set of (68) constraints, for example an
eigenvector expansion of the matrix $\sum_{k=0}^{n-1}u_{jk}u_{j^{\prime}k}$
followed by eigenvalues adjustment[17], Gram–Schmidt orthogonalization, etc.
However, the SVD expansion (74) is special, see (118) below.

  3. 3.

Put this new $u_{jk}$ to (71a), then multiply it by $u^{*}_{ik}$ and sum over
$k=0\dots n-1$. As the $u_{jk}$ is partially unitary (68) obtain new values
for Lagrange multipliers $\widetilde{\lambda}_{ij}$ and take it’s Hermitian
part999The equation for Lagrange multipliers (76) produces an arbitrary matrix
$\widetilde{\lambda}_{ij}$; a variation of the constraints produces Hermitian
matrix $\lambda_{ij}$. Lagrange multipliers in (71) should be set to make the
first variation at given $u_{jk}$ as close to zero as possible; least squares
expansion of the first variation ($D\times n$ matrix) in Lagrange multipliers
($D\times D$ matrix) gives (76). For an arbitrary matrix $\mathcal{A}$ it’s
best approximation by a Hermitian matrix $\mathcal{B}$ is the Hermitian part
$\mathcal{B}=\mathcal{A}_{H}=\frac{1}{2}(\mathcal{A}+\mathcal{A}^{\dagger})$.
This follows immediately from the Frobenius norm triangle inequality by
splitting the matrix into Hermitian $\mathcal{A}_{H}$ and anti–Hermitian
$\mathcal{A}_{AH}$ parts:
$\|\mathcal{B}-\mathcal{A}\|_{F}=\|\mathcal{B}-\mathcal{A}_{H}-\mathcal{A}_{AH}\|_{F}\leq\|\mathcal{B}-\mathcal{A}_{H}\|_{F}+\|\mathcal{A}_{AH}\|_{F}$.
:

$\displaystyle\widetilde{\lambda}_{ij}$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{j^{\prime}k^{\prime}}S_{j^{\prime}k^{\prime};jk}u^{*}_{ik}$
(76) $\displaystyle\lambda_{ij}$
$\displaystyle=\frac{1}{2}\left[\widetilde{\lambda}_{ij}+\widetilde{\lambda}^{*}_{ji}\right]$
$\displaystyle i,j=0\dots D-1$ (77)

This $\lambda_{ij}$ is the next iteration of Lagrange multipliers. As
iterations proceed – the $\widetilde{\lambda}_{ij}$ is expected to converge to
a Hermitian matrix by itself, without (77) required. For original (not yet
full–constraint adjusted) $u_{jk}$, which is an eigenvector of
$S_{j^{\prime}k^{\prime};jk}$, the $\widetilde{\lambda}_{ij}$ is Hermitian.
The anti–Hermitian part of $\widetilde{\lambda}_{ij}$ cancels in the quadratic
form (70). One can possibly obtain a Hermitian $\lambda_{ij}$ right away with
multiplication of (71a) by itself (instead of $u_{jk}$ for (76)); the
Hermitian $\lambda_{ij}$ is then obtained from $\lambda_{ij}^{2}$ as all the
eigenvalues of $\lambda_{ij}$ are all positive; the result is very similar to
(77), a drawback for this new $\lambda_{ij}$ — the (73) now holds only
approximately for current iteration of $u_{jk}$, see
com/polytechnik/utils/KGOIterationalLambda2.java.

  4. 4.

Put this new $\lambda_{ij}$ to (70) and repeat iteration process until
converged. On the first iteration take initial values of Lagrange multipliers
$\lambda_{ij}=0$.

For a simpler scalar QCQP optimization problem of [15], ‘‘ _Appendix F:
Directional Information: $I\xrightarrow[\psi]{\quad}\max$ Subject To the
Constraint $\Braket{\psi}{C}{\psi}=0$_’’, where we considered a single
quadratic constraint, similar iteration algorithm converges fast but may fail
when optimization and constraint matrices have a number of eigenvectors in
common. The optimization problem (67) subject to (68) constraints is a problem
of (52) type, it has a more complex internal structure than the problem
considered in [15].

The described Lagrange multipliers algorithm is based on eigenvalue problem
solution: (70) with partial constraint (69) as normalizing:
$D=\Braket{\psi^{2}}$. It is much less sensitive to degeneracy than Newtonian
type iterations, where even a single degenerate degree of freedom makes linear
system (with Hessian matrix) iteration to fail. A question arise when the
described above iteration algorithm fails. Currently — we do not have the
exact answer; the condition of iteration algorithm convergence requires a
separate study. The algorithm does not converge well for partially unitary
operators with $D<n$, but given large enough iterations number it produces a
good enough solution. The reason for a slow convergence is that with (76)
$\lambda_{ij}$ the Hessian matrix is degenerated at the adjusted $u_{jk}$ (75)
— at this $u_{jk}$ not only first but also second variation of the objective
function is zero; this is a constraint qualification problem. The algorithm
does not diverge, it provides a sequence of close to optimal solutions. See
com/polytechnik/utils/KGOIterationalSimpleOptimizationU.java for a numerical
implementation. We also tried to find an algorithm of contraction mapping
type, but this requires more study. The convergence can be greatly improved
using linear constraints, see Appendix A.4 below where the constraints (68)
were replaces by the closeness of $u_{jk}$ to current iteration value (106).
In many situation, however, an approximate solution is sufficient.

###  A.1 On Constrained Optimization In The Singular Values Basis

Before we go further let us discuss the roles of (74) singular values and
their relation to the calculation of Lagrange multipliers. If we write
optimization problem (70) in SVD basis (74) the $u_{jk}$ is represented as a
product of three matrices. The constraints (68) require all singular values
$\Sigma_{jj}=\pm 1$. We denote this diagonal matrix as vector $\Sigma_{s}$.
The objective function (65) is then
$\mathcal{F}=\sum_{s,s^{\prime}=0}^{D-1}\Sigma_{s}\widetilde{S}_{ss^{\prime}}\Sigma_{s^{\prime}}$.
Obtain constrained optimization problem with $D$ Lagrange multipliers
$\widetilde{\lambda}_{s}$

$\displaystyle\sum\limits_{s,s^{\prime}=0}^{D-1}\Sigma_{s}\widetilde{S}_{ss^{\prime}}\Sigma_{s^{\prime}}+\sum\limits_{s=0}^{D-1}\widetilde{\lambda}_{s}\left[1-\Sigma_{s}^{2}\right]\to\max$
(78)
$\displaystyle\widetilde{S}_{ss^{\prime}}=\sum\limits_{j,j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}U_{js}V^{\dagger}_{sk}S_{jk;j^{\prime}k^{\prime}}U_{j^{\prime}s^{\prime}}V^{\dagger}_{s^{\prime}k^{\prime}}$
(79)

from which we immediately obtain the values

$\displaystyle\widetilde{\lambda}_{s}$
$\displaystyle=\frac{1}{\Sigma_{s}}\sum\limits_{s^{\prime}=0}^{D-1}\widetilde{S}_{ss^{\prime}}\Sigma_{s^{\prime}}$
(80)

for all adjusted $\Sigma_{s}=1$

$\displaystyle\widetilde{\lambda}_{s}$
$\displaystyle=\sum\limits_{s^{\prime}=0}^{D-1}\widetilde{S}_{ss^{\prime}}$
(81)

Comparing (78) with (70) obtain $\lambda_{ij}$ in original basis

$\displaystyle\widetilde{\lambda}_{s}$
$\displaystyle=\sum\limits_{i,j=0}^{D-1}\lambda_{ij}U_{is}U_{js}$ (82)
$\displaystyle\lambda_{ij}$
$\displaystyle=\sum\limits_{s=0}^{D-1}U_{is}U_{js}\widetilde{\lambda}_{s}$
(83)

Whereas the original functional (70) has $D^{2}$ Lagrange multipliers
$\lambda_{ij}$, the (78) has only $D$ — a constraint for every singular value
of the matrix $u_{jk}$; it is clear why: since the partial constraint (69) is
always satisfies from the eigenproblem it is sufficient to set $D-1$ diagonal
elements of (68) to $1$, then all off–diagonal elements are immediately zero.

###  A.2 On Iteration Step Without Using The SVD

In the algorithm above we extensively used SVD expansion (74) for iterations.
Let us consider how to avoid using the SVD by replacing it with an eigenvalue
problem of the dimension $D\times D$ for the purpose of both: computational
complexity and better understanding of the algorithm. Obtained partial
constraint (69) solution matrix $u_{jk}$ is non–orthogonal, the Gram matrix
is:

$\displaystyle G^{u}_{jj^{\prime}}$
$\displaystyle=\sum\limits_{k=0}^{n-1}u_{jk}u_{j^{\prime}k}$ (84)

We need to ‘‘adjust’’ $u_{jk}$ to satisfy the full set of (68) constraints.
Consider the eigenstates of the Gram matrix

$\displaystyle\Ket{G^{u}|u^{[i]}}$
$\displaystyle=\lambda_{G}^{[i]}\Ket{u^{[i]}}$ (85)

The eigenvalues of this problem are equal to the singular values (74) squared
$\lambda_{G}^{[i]}=\Sigma^{2}_{ii}$. Setting all $\lambda_{G}^{[i]}=1$
(eigenvalues adjustment technique [17]) produces a new basis in which (68)
constraints are satisfied in full. The result is identical to the transform
(75) of setting all $\Sigma_{jj}=1$ but it is obtained without solving a SVD
problem, the eigenvalue $D\times D$ problem (85) is used instead, see
com/polytechnik/utils/KGOEVSelection.java:getEVAdjustedTo1() for an
implementation.

Optimization problem is question is invariant relatively a unitary transform
(the $A_{sj}$ is a unitary matrix)

$\displaystyle v_{sk}$ $\displaystyle=\sum\limits_{j=0}^{D-1}A_{sj}u_{jk}$
(86)

The tensor $S_{jk;j^{\prime}k^{\prime}}$ transforms with $A_{sj}$ as (87),
Gram matrix (84) corresponds to the tensor
$S_{jk;j^{\prime}k^{\prime}}=G^{u}_{jj^{\prime}}\delta_{kk^{\prime}}$.

$\displaystyle S_{sk;s^{\prime}k^{\prime}}$
$\displaystyle=\sum\limits_{jj^{\prime}=0}^{D-1}A_{sj}S_{jk;j^{\prime}k^{\prime}}A_{s^{\prime}j^{\prime}}$
(87) $\displaystyle\mathcal{F}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}S_{jk;j^{\prime}k^{\prime}}u_{j^{\prime}k^{\prime}}=\sum\limits_{s,s^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}v_{sk}S_{sk;s^{\prime}k^{\prime}}v_{s^{\prime}k^{\prime}}$
(88)

The constraints for new variables $v_{sk}$ have the same form (68)

$\displaystyle\delta_{ss^{\prime}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}v_{sk}v_{s^{\prime}k}$
$\displaystyle s,s^{\prime}=0\dots D-1$ (89)

Let us transform the input to the basis of Gram matrix eigenvectors. Solve
generalized eigenproblem (85) to find the eigenvalues $\lambda_{G}^{[s]}$ and
the eigenvectors $v^{[s]}_{j}$ of the Gram matrix $G^{u}_{jj^{\prime}}$

$\displaystyle\sum\limits_{j^{\prime}=0}^{D-1}G^{u}_{jj^{\prime}}v^{[s]}_{j^{\prime}}$
$\displaystyle=\lambda_{G}^{[s]}v^{[s]}_{j}$ (90)

Were it all $\lambda_{G}^{[s]}=1$ — the eigenstates of the Gram matrix would
form the sought partially unitary operator, but this is typically not. Take
Gram matrix eigenvectors as a new basis, the unitary transform matrix is
$A_{sj}=v^{[s]}_{j}$, and write optimization problem (88) in this new basis
$v_{sk}$ (86) with the tensor $S_{sk;s^{\prime}k^{\prime}}$ transformed from
the $S_{jk;j^{\prime}k^{\prime}}$ according to (87). If all scaling
coefficients $\mu_{s}=1$ — this would be exactly the original problem since it
is invariant relatively unitary transforms of the basis, but if we put the
factors $\mu_{s}$ (91) — this makes the solution to satisfy (89); non–unitary
scaling factors $\mu_{s}$ adjust the solution to satisfy the full set of the
constraints.

$\displaystyle\mu_{s}=\pm\frac{1}{\sqrt{\lambda_{G}^{[s]}}}$ (91)
$\displaystyle\sum\limits_{s,s^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}\mu_{s}v_{sk}S_{sk;s^{\prime}k^{\prime}}v_{s^{\prime}k^{\prime}}\mu_{s^{\prime}}\xrightarrow[v]{\quad}\max$
(92)

This scaling adjustment performed in Gram matrix basis is an alternative to
SVD adjustment (75). One need to convert the problem from original basis to
the basis of Gram matrix eigenvectors, then scale them by the (91) factors.
The $\mu_{s}v_{sk}$ satisfies partial orthogonality constraints. We can write
optimization problem in this new basis, and perform the iteration algorithm of
Appendix A above, then ‘‘chaining’’ unitary transforms as iterations proceed,
the result will be identical as the problem is invariant relatively these
transforms, but the idea of solution adjustment in the from of pure scaling
opens a number of new ways to improve the algorithm, see
com/polytechnik/utils/KGOIterationalMultipleTransforms.java for a numerical
implementation.

###  A.3  On Operator–Dependent Solution Adjustment

In the previous section we considered solution adjustment procedure applied to
some initial ‘‘partial’’ solution. This adjustment is a non–unitary basis
transform. A question arise about a generalization: applying some other
non–unitary transform before the adjustment. Optimization problem in question
is to maximize (93) subject to (94) constraints:

$\displaystyle\frac{\mathcal{F}}{D}=$
$\displaystyle\frac{\sum\limits_{j,j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}S_{jk;j^{\prime}k^{\prime}}u_{j^{\prime}k^{\prime}}}{\sum\limits_{j=0}^{D-1}\sum\limits_{k=0}^{n-1}u^{2}_{jk}}\xrightarrow[u]{\quad}\max$
(93) $\displaystyle\delta_{jj^{\prime}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}u_{j^{\prime}k}$
$\displaystyle j,j^{\prime}=0\dots D-1$ (94)

Consider a Hermitian operator $\mathcal{J}$ with matrix elements
$\mathcal{J}_{jj^{\prime}}$, this can be e.g. Lagrange multipliers matrix
(77), unit matrix, etc. A generalized eigenvalue problem with
$\mathcal{J}_{jj^{\prime}}$ and $G^{u}_{jj^{\prime}}$ (84) matrices is
formulated as

$\displaystyle\sum\limits_{j^{\prime}=0}^{D-1}\mathcal{J}_{jj^{\prime}}v^{[s]}_{j^{\prime}}$
$\displaystyle=\lambda_{\mathcal{J}}^{[s]}\sum\limits_{j^{\prime}=0}^{D-1}G^{u}_{jj^{\prime}}v^{[s]}_{j^{\prime}}$
(95)

Because of the Gram matrix $G^{u}_{jj^{\prime}}$ in the right hand side
obtained solution

$\displaystyle v_{sk}$
$\displaystyle=\sum\limits_{j=0}^{D-1}v^{[s]}_{j}u_{jk}$ (96)

satisfies (94) constraints
$\delta_{ss^{\prime}}=\sum_{k,k^{\prime}=0}^{n-1}v_{sk}v_{s^{\prime}k}$. The
transform $v^{[s]}_{j}$ is non–unitary

$\displaystyle\delta_{ss^{\prime}}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{D-1}v^{[s]}_{j}G^{u}_{jj^{\prime}}v^{[s^{\prime}]}_{j^{\prime}}$
(97a) $\displaystyle G^{u;-1}_{jj^{\prime}}$
$\displaystyle=\sum\limits_{s=0}^{D-1}v^{[s]}_{j}v^{[s]}_{j^{\prime}}$ (97b)

Condition (97a) creates the basis (96) satisfying partial orthogonality
constraints. Let us write the optimization problem (93) in this new basis
$v_{sk}$. Using

$\displaystyle u_{jk}$
$\displaystyle=\sum\limits_{j^{\prime},s=0}^{D-1}G^{u}_{jj^{\prime}}v^{[s]}_{j^{\prime}}v_{sk}$
(98)

obtain the original problem (93) with the tensor $S_{sk;s^{\prime}k^{\prime}}$
instead of $S_{jk;j^{\prime}k^{\prime}}$

$\displaystyle S_{sk;s^{\prime}k^{\prime}}$
$\displaystyle=\sum\limits_{j,j^{\prime},i,i^{\prime}=0}^{D-1}v^{[s]}_{i}G^{u}_{ij}S_{jk;j^{\prime}k^{\prime}}G^{u}_{i^{\prime}j^{\prime}}v^{[s^{\prime}]}_{i^{\prime}}$
(99)

This is a generalization of (87) to non–unitary transforms. This is exactly
the original problem (without an adjustment), but written in the $v_{sk}$
basis.

It can be noticed that adjustment procedure of previous section is actually a
non–unitary transform with the inverse square root of the Gram matrix
$G^{u;-1/2}_{jj^{\prime}}$ (91); there are $2^{D-1}$ distinct combinations of
signs but we take all equal to $1$. The adjustment is equivalent to
multiplying (98) by $G^{u;-1/2}_{jj^{\prime}}$ to obtain the ‘‘adjusted’’
tensor

$\displaystyle S^{adj}_{sk;s^{\prime}k^{\prime}}$
$\displaystyle=\sum\limits_{j,j^{\prime},i,i^{\prime}=0}^{D-1}v^{[s]}_{i}G^{u;1/2}_{ij}S_{jk;j^{\prime}k^{\prime}}G^{u;1/2}_{i^{\prime}j^{\prime}}v^{[s^{\prime}]}_{i^{\prime}}$
(100)

This way the adjustment is ‘‘transferred’’ from the state $u_{jk}$ to operator
$S_{sk;s^{\prime}k^{\prime}}$. Equivalent $u_{jk}$ adjustment corresponds to
$u_{jk}=\sum_{j^{\prime},s=0}^{D-1}G^{u;1/2}_{jj^{\prime}}v^{[s]}_{j^{\prime}}v_{sk}$.
The (100) is an important option to transfer an adjustment from a state to
tensor, this allows to combine the adjustment with optimization algorithm.
Considered in Section A.2 above adjustment procedure corresponds to
$\mathcal{J}_{jj^{\prime}}$ being a unit matrix.

###  A.4  On Optimization Algorithm With Linear Constraint Iteration

In previous sections we considered optimization algorithm with quadratic
constraints of (68) form. In numerical implementation (70) these constraints
lead to a poor convergence since at the point the constraints are applied the
Hessian matrix is degenerated. Consider a linear type of constraints.

Extend $u_{jk}$ with one more degree of freedom $\chi$ to form a ‘‘vector’’ of
the dimension $D\times n+1$.

$\displaystyle\bm{z}$ $\displaystyle=\begin{pmatrix}u_{jk}\\\
\chi\end{pmatrix}$ (101)

Then the quadratic from

$\displaystyle\mathcal{F}$
$\displaystyle=\frac{\bm{z}^{T}\mathcal{S}\bm{z}}{\bm{z}^{T}\bm{z}}=\frac{\chi^{2}S_{0}+2\chi\sum\limits_{j=0}^{D-1}\sum\limits_{k=0}^{n-1}b_{jk}u_{jk}+\sum\limits_{j,j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}S_{jk;j^{\prime}k^{\prime}}u_{j^{\prime}k^{\prime}}}{\chi^{2}+\sum\limits_{j=0}^{D-1}\sum\limits_{k=0}^{n-1}u^{2}_{jk}}$
(102)

has the matrix $\mathcal{S}$

$\displaystyle\mathcal{S}$
$\displaystyle=\left(\begin{array}[]{c|c}S_{jk;j^{\prime}k^{\prime}}&b_{j^{\prime}k^{\prime}}\\\
\hline\cr b_{jk}&S_{0}\end{array}\right)$ (105)

The idea is to consider the $b_{jk}$ and $S_{0}$ as some kind of ‘‘Lagrange
Multipliers’’ to set the variation of (102) to zero at the ‘‘adjusted’’
$u_{jk}$, denote it as iteration value $u^{IT}_{jk}$. Consider the constraints

$\displaystyle u_{jk}-u^{IT}_{jk}=0$ (106)

this is the closeness of $u_{jk}$ to current iteration value $u^{IT}_{jk}$
(adjusted value satisfying all the required constraints (68)). A one more
degree of freedom $\chi$ was introduced to preserve the form of the Rayleigh
quotient for the optimization problem (102). Variating it over $u_{jk}$ and
$\chi$ obtain (109) and (110) respectively; in these formulas
$u_{jk}=u^{IT}_{jk}$, $F_{0}^{IT}$ is a known constant, and $B_{0}$ and
$S_{0}$ are unknown constants.

$\displaystyle F_{0}^{IT}$
$\displaystyle=\sum\limits_{j,j^{\prime}=0}^{D-1}\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}S_{jk;j^{\prime}k^{\prime}}u_{j^{\prime}k^{\prime}}$
(107) $\displaystyle B_{0}$
$\displaystyle=\sum\limits_{j=0}^{D-1}\sum\limits_{k=0}^{n-1}b_{jk}u_{jk}$
(108) $\displaystyle 0$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{D-1}\sum\limits_{k^{\prime}=0}^{n-1}S_{jk;j^{\prime}k^{\prime}}u_{j^{\prime}k^{\prime}}+\chi
b_{jk}-\frac{\chi^{2}S_{0}+2\chi B_{0}+F_{0}^{IT}}{D+\chi^{2}}u_{jk}$ (109)
$\displaystyle 0$ $\displaystyle=\chi S_{0}+B_{0}-\frac{\chi^{2}S_{0}+2\chi
B_{0}+F_{0}^{IT}}{D+\chi^{2}}\chi$ (110) $\displaystyle 0$
$\displaystyle=F_{0}^{IT}+\chi B_{0}-\left(\chi^{2}S_{0}+2\chi
B_{0}+F_{0}^{IT}\right)\frac{D}{D+\chi^{2}}$ (111)

Multiply (109) by $u_{jk}$ and sum it over $j$ and $k$, obtain (111). For a
given $\chi$ the (110) and (111) can be considered as a linear system for
$B_{0}$ and $S_{0}$. Obtained $2\times 2$ linear system is degenerated and has
multiple solutions:

$\displaystyle(D-\chi^{2})B_{0}+(D\chi)S_{0}$ $\displaystyle=F_{0}^{IT}\chi$
(112)

The specific set $(\chi,S_{0},B_{0})$ should be selected for best convergence.
The selection

$\displaystyle\chi$ $\displaystyle=1$ (113a) $\displaystyle S_{0}$
$\displaystyle=F_{0}^{IT}$ (113b) $\displaystyle B_{0}$ $\displaystyle=-S_{0}$
(113c)

is the first one to try.

  1. 1.

Take the values of $b_{jk}$ and $S_{0}$ to construct (105).

  2. 2.

Solve (102) and select the most appropriate vector $\mathbf{z}$. The result of
this step — the ‘‘adjusted’’ $u^{IT}_{jk}$ satisfying all the required
constraints (68).

  3. 3.

Take this new $u^{IT}_{jk}$, and select some value of $\chi$, for example
(113), calculate ‘‘Lagrange Multipliers’’ $b_{jk}$ (109) and $S_{0}$ (112) to
construct (105) matrix. If one uses $\chi$ value from (102) maximization
problem — iterations typically stick to some local maximum. If one uses a
fixed value for $\chi$, such as (113a) — a convergence is observed; not very
fast, but better than in the Appendix A above. Repeat iteration process. On
the first iteration take $b_{jk}=S_{0}=0$.

This ‘‘Linear constraints’’ algorithm is implemented in the
com/polytechnik/utils/KGOIterationalLinearConstraintsE.java. An attempt to use
$D$ extra degrees of freedom instead of a single one was much less successful
com/polytechnik/utils/KGOIterationalLinearConstraintsExtraDegreesOfFreedom.java.

###  A.5  An Algorithm to Find an Approximate Solution to the Knowledge
Generalizing Operator

Consider the same problem (67) subject to (68) constraint. The bases are
considered already orthogonalized:
$\delta_{jj^{\prime}}=\Braket{f_{j}}{f_{j^{\prime}}}$ and
$\delta_{kk^{\prime}}=\Braket{x_{k}}{x_{k^{\prime}}}$. Assume we found
optimization problem (67) solution with ‘‘partial’’ constraints (69), this is
(70) with $\lambda_{ij}=0$. Put it to (66) and expand $u_{jk}$ in SVD:

$\displaystyle u_{jk}$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{D-1}\sum\limits_{k^{\prime}=0}^{n-1}U_{jj^{\prime}}\Sigma_{j^{\prime}k^{\prime}}V^{\dagger}_{k^{\prime}k}$
(114) $\displaystyle\Braket{f_{j}}{f_{j^{\prime}}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}u_{jk}\Braket{x_{k}}{x_{k^{\prime}}}u^{*}_{j^{\prime}k^{\prime}}$
$\displaystyle j,j^{\prime}=0\dots D-1$ (115)

Write (115) for orthogonal bases $\mathfrak{x}_{k}$ and $\mathfrak{f}_{j}$

$\displaystyle\mathfrak{x}_{k}$
$\displaystyle=\sum\limits_{k^{\prime}=0}^{n-1}V^{\dagger}_{kk^{\prime}}x_{k^{\prime}}$
(116) $\displaystyle\mathfrak{f}_{j}$
$\displaystyle=\sum\limits_{j^{\prime}=0}^{D-1}U^{\dagger}_{jj^{\prime}}f_{j^{\prime}}$
(117) $\displaystyle\Braket{\mathfrak{f}_{j}}{\mathfrak{f}_{j^{\prime}}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}\Sigma_{jk}\Braket{\mathfrak{x}_{k}}{\mathfrak{x}_{k^{\prime}}}\Sigma_{j^{\prime}k^{\prime}}$
$\displaystyle j,j^{\prime}=0\dots D-1$ (118)

The (118) is (115) written in $\mathfrak{x}_{k}$ and $\mathfrak{f}_{j}$
orthogonal bases. Since
$\delta_{jj^{\prime}}=\Braket{\mathfrak{f}_{j}}{\mathfrak{f}_{j^{\prime}}}$
and
$\delta_{kk^{\prime}}=\Braket{\mathfrak{x}_{k}}{\mathfrak{x}_{k^{\prime}}}$
the (118) is satisfied only when all singular values of $u_{jk}$ are $\pm 1$.
Actually we made a single iteration of the algorithm above, this
$\Sigma_{jk}$-adjusted solution is an approximate solution one should try
first. Since $\Sigma_{jk}$ is diagonal, in the
$(\mathfrak{f}_{j},\mathfrak{x}_{k})$ basis we have a one–to–one relation

$\displaystyle\mathfrak{f}_{j}$ $\displaystyle=\mathfrak{x}_{j}\Sigma_{jj}$
(119)

This is not a least squares type of relation, for example the result is
invariant relatively the transform $\Sigma_{jk}\to-\Sigma_{jk}$. If
$(\mathfrak{f}_{j},\mathfrak{x}_{k})$ basis satisfies (118) then all singular
values are $\pm 1$ (the condition of partial unitarity) and
$(\mathfrak{f}_{j},\mathfrak{x}_{k})$ relation is plain
$\mathfrak{f}_{j}=\pm\,\mathfrak{x}_{j}$. The probability (37) in this basis
is

$\displaystyle\Braket{\psi_{\mathfrak{f}}}{\mathcal{U}}{\psi_{\mathfrak{x}}}^{2}$
$\displaystyle=\frac{\left|\sum\limits_{j=0}^{D-1}\mathfrak{f}_{j}\mathfrak{x}_{j}\Sigma_{jj}\right|^{2}}{\sum\limits_{j=0}^{D-1}\mathfrak{f}^{2}_{j}\sum\limits_{k=0}^{n-1}\mathfrak{x}^{2}_{k}}$
(120)

Partial unitarity ‘‘adjusted’’ case corresponds to $\Sigma_{jj}=\pm 1$.

Consider the meaning of a state with an arbitrary $\Sigma_{jj^{\prime}}$. The
(118) is actually the constraint (68) but with the positive diagonal matrix
$\Sigma^{2}_{jj^{\prime}}$, not $\delta_{jj^{\prime}}$. What does it mean if
we put this $u_{jk}$ ‘‘partial constraint (69)’’ solution to probability (37)
without any adjustment? This breaks the preservation of probability, the
probability (120) is no longer $[0:1]$ bounded, it is now $0\leq
P(\mathbf{f})\Big{|}_{\mathbf{x}}\leq\max\limits_{j}\Sigma^{2}_{jj}$; the
range $[0:1]$ holds only ‘‘on average’’, for the entire sample. However, this
does not change the calculation of outcome value (49). One can also modify
(120) to have the probability $[0:1]$ bounded, the maximal value is $1$, it
corresponds to $\mathfrak{f}_{j}=\mathfrak{x}_{j}\Sigma_{jj}$.

$\displaystyle P(\mathfrak{f})\Big{|}_{\mathfrak{x}}$
$\displaystyle\approx\frac{\left|\sum\limits_{j=0}^{D-1}\mathfrak{f}_{j}\mathfrak{x}_{j}\Sigma_{jj}\right|^{2}}{\sum\limits_{j=0}^{D-1}\mathfrak{f}^{2}_{j}\sum\limits_{j=0}^{D-1}\mathfrak{x}^{2}_{j}\Sigma^{2}_{jj}}$
(121)

But this is only for evaluation, this is not the function used in optimization
problem, optimization problem with the probability (122) is much more
difficult. There is a trivial option to use the probability (123) for
optimization and (121) for evaluation. The $\Sigma^{2}_{jj}$, $j=0\dots D-1$,
factor (whether the singular values are adjusted or not) in the denominator
prevents a decrease of probability when polluting the $\mathbf{x}$–space with
a large number of completely random components (127); the value of
$\mathfrak{f}$ (49) does not depend on this $\mathfrak{x}$–depended factor,
maximal value of probability corresponds to
$\mathfrak{f}_{j}=\mathfrak{x}_{j}\Sigma_{jj}$; the probability is invariant
with respect to $\mathfrak{f}_{j}\to C\mathfrak{f}_{j}$, normalize it to const
to obtain actual values. This partial constraint solution of (67) subject to
(69) is an approximate solution one may try. Whereas a quantum channel that
preserves probability ‘‘on average’’ does not have a physical meaning, in data
analysis it is an approximation with a clear meaning: it emphasizes (118)
internal relations with high probability, the $\Sigma_{jj}$ factor in (121).
Mathematically this means that in (52) we allow operators $\mathcal{U}$ that
preserve Gram matrix spur: $D=\sum\limits_{j=0}^{D-1}\Sigma^{2}_{jj}$, not the
Gram matrix itself (115) as previously considered; the solution can be found
from eigenproblem (70) in original basis with $\lambda_{ij}=0$.

Conceptually, this technique consists in taking any approximate $u_{jk}$, such
as least squares (7) or any other matrix, not necessary (70) solution, Gram
matrix spur preservation is not required, expanding $u_{jk}$ in SVD (114),
then set $\Sigma_{jj}$ to $1$ or $-1$. There are $2^{m-1}$ distinct
combinations, typically the minimal change adjustment — all $\Sigma_{jj}=1$
gives the best result as the initial $\Sigma_{jj}$ are positive. Obtained new
$u_{jk}$ matrix with singular values equal to $\pm 1$ satisfies all the
required constraints (115). Alternatively one can solve the eigenproblem (85)
and adjust all the eigenvectors by the factors $\pm 1/\sqrt{\lambda^{[j]}}$
(91) to obtain the same solution without using the SVD (it is equivalent to
multiplication (100) of unadjusted $u_{jk}$ by inverse square root of
corresponding Gram matrix).

##  Appendix B On Adjusted Normalizing Of Probability

The probability (37) has a normalizing factor as a product of two Christoffel
functions: on $\mathbf{x}$ and on $\mathbf{f}$ (33); these two Christoffel
functions have $n$ and $m$ degrees of freedom respectively. In some situations
it is convenient to construct a normalizing factor where both $\mathbf{x}$\-
and $\mathbf{f}$\- factors have the same number of degrees of freedom: $m$.

One can consider the probability adjusted to only ‘‘important’’
$\mathbf{x}$–components, this is
$\Braket{\psi_{\mathbf{g}}}{\psi_{u(\mathbf{y})}}^{2}$ from (64) expanded:

$\displaystyle\Braket{\psi_{\mathbf{g}}}{\mathcal{U}}{\psi_{\mathbf{y}}}^{2}$
$\displaystyle=\frac{\left|\sum\limits_{k=0}^{n-1}\sum\limits_{j,s=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{js}u_{sk}y_{k}\right|^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}\sum\limits_{j,j^{\prime}=0}^{m-1}\sum\limits_{k,k^{\prime}=0}^{n-1}y_{k}u_{jk}G^{\mathbf{f};\,-1}_{jj^{\prime}}u_{j^{\prime}k^{\prime}}y_{k^{\prime}}}$
(122)

Whereas this formula for
$\Braket{\psi_{\mathbf{g}}}{\psi_{u(\mathbf{y})}}^{2}$ has a more suitable
normalizing than (37), it has $u_{jk}$ in the denominator and the problem can
no longer101010 This difficulty does not arise with $\mathbf{x}$\- and
$\mathbf{f}$\- being the same space. For example for a unitary $\mathcal{U}$
the denominator does not depend on $\mathcal{U}$.  be reduced to the one of
form (40) that requires only the moments of Christoffel functions product
(33). For probability evaluation, not for optimization, this can be done
straightforward (121). A quantum channel $u_{jk}$ optimizing (18) with the
probability (122) is an interesting direction of future research, this new
problem is no longer a QCQP problem — it is a problem to maximize the sum of
$M$ ratios of two quadratic forms on $u_{jk}$ subject to (38) constraint or,
more generally, an unconstrained optimization of (64). The one in the
nominator is a dyadic product squared, the one in the denominator is
non–negative, it cancels with the nominator when it’s value is close to zero.

To adjust the number of degrees of freedom one can use a much simpler
alternative approach. All we need is to calculate a Christoffel function in
$\mathbf{x}$–space to normalize the probability. A trivial approach is to use
the contributing subspace $\Ket{\phi^{[i]}}$, e.g. from (20). Despite the
moments $\Braket{f_{j}x_{k}}$ have been used to build the contributing
subspace $\Ket{\phi^{[i]}}$, this does not create any difficulty as we use
these projections only to construct a Christoffel function with matched number
of degrees of freedom. The (126) is invariant with respect to
$\Braket{f_{j}x_{k}}\to-\Braket{f_{j}x_{k}}$ and tends to a constant when any
$\Braket{f_{j}x_{k}}\to\infty$ (factors in the denominator and inverse
$G^{fxa}_{jj^{\prime}}$ matrix).

$\displaystyle\Braket{\psi_{\mathbf{g}}}{\mathcal{U}}{\psi_{\mathbf{y}}}^{2}$
$\displaystyle=\frac{\left|\sum\limits_{k=0}^{n-1}\sum\limits_{j,s=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{js}u_{sk}y_{k}\right|^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}\sum\limits_{i=0}^{m-1}{\phi^{[i]}}^{2}(\mathbf{y})}$
(123) $\displaystyle G^{fxa}_{jj^{\prime}}$
$\displaystyle=\sum\limits_{k,k^{\prime}=0}^{n-1}\Braket{f_{j}x_{k}}G^{\mathbf{x};\,-1}_{kk^{\prime}}\Braket{x_{k^{\prime}}f_{j^{\prime}}}$
(124) $\displaystyle G^{C}_{qq^{\prime}}$
$\displaystyle=\sum\limits_{s,s^{\prime}=0}^{n-1}\sum\limits_{j,j^{\prime}=0}^{m-1}G^{\mathbf{x};\,-1}_{qs}\Braket{x_{s}f_{j}}G^{fxa;\,-1}_{jj^{\prime}}\Braket{f_{j^{\prime}}x_{s^{\prime}}}G^{\mathbf{x};\,-1}_{s^{\prime}q^{\prime}}$
(125) $\displaystyle K^{adj}(\mathbf{x})$
$\displaystyle=\frac{1}{\sum\limits_{i=0}^{m-1}{\phi^{[i]}}^{2}(\mathbf{x})}=\frac{1}{\sum\limits_{q,q^{\prime}=0}^{n-1}x_{q}G^{C}_{qq^{\prime}}x_{q^{\prime}}}$
(126)

The value of $K^{adj}(\mathbf{x})$ is never zero on training sample since
contributing subspace always has a constant among the components. The
probability (123) uses Christoffel function with adjusted number of degrees of
freedom $K^{adj}(\mathbf{x})$ (126) instead of the original $K(\mathbf{x})$
(10) for the probability (37). The difference between two these Christoffel
functions is in extra terms in the denominator sum. Since the entire
$\mathbf{x}$–space can be represented as the direct sum of $\Ket{\phi^{[i]}}$
and $\Ket{\phi^{O;[i]}}$, a subspace of $\mathbf{x}$ orthogonal to
$\Ket{\phi^{[i]}}$, the $K(\mathbf{x})$ (10) is:

$\displaystyle K(\mathbf{x})$
$\displaystyle=\frac{1}{\sum\limits_{i=0}^{m-1}{\phi^{[i]}}^{2}(\mathbf{x})+\sum\limits_{i=m}^{n-1}{\phi^{O;[i]}}^{2}(\mathbf{x})}=\frac{1}{\sum\limits_{k,k^{\prime}=0}^{n-1}x_{k}G^{\mathbf{x};\,-1}_{kk^{\prime}}x_{k^{\prime}}}$
(127)

Thus we always have $K^{adj}(\mathbf{x})\geq K(\mathbf{x})$. The moments of
two Christoffel functions product are

$\displaystyle\Braket{x_{k}f_{j}}{K^{adj(\mathbf{x})}K^{(\mathbf{f})}}{x_{k^{\prime}}f_{j^{\prime}}}$
$\displaystyle=\sum\limits_{l=0}^{M}\omega^{(l)}\frac{x^{(l)}_{k}x^{(l)}_{k^{\prime}}}{\sum\limits_{q,q^{\prime}=0}^{n-1}x_{q}^{(l)}G^{C}_{qq^{\prime}}x_{q^{\prime}}^{(l)}}\cdot\frac{f^{(l)}_{j}f^{(l)}_{j^{\prime}}}{\sum\limits_{s,s^{\prime}=0}^{m-1}f^{(l)}_{s}G^{\mathbf{f};\,-1}_{ss^{\prime}}f^{(l)}_{s^{\prime}}}$
(128)

This tensor has the same dimensions as (33), the difference only in
normalizing — it uses $G^{C}_{qq^{\prime}}$ from (125) instead of
$G^{\mathbf{x};\,-1}_{qq^{\prime}}$ in (33). Despite it now depends on
$\Braket{f_{j}x_{k}}$ moments — they are used only to construct Christoffel
function for normalizing, this does not change the essence of the solution due
to the invariance properties of the Christoffel function.

##  Appendix C On Contributing Subspace Selection

Considered above optimization problem finds partially unitary operator
$u_{jk}$ that does both: selects the contributing subspace ($m$ vectors of the
dimension $n$) and optimizes the objective function. Besides computational
difficulties this also creates a problem with normalizing since properly
normalized objective function (64) has operator $u_{jk}$ both in the nominator
and in the denominator (122), thus some surrogate normalizing (39), (128),
(46), or (47) was used instead. It is a very attractive option to split the
problem into two:

  * •

Find the contributing subspace $\phi^{[j]}$ of the dimension $m$.

  * •

Find a unitary (not partially unitary!) operator $\mathcal{U}$ mapping from
$\phi^{[j]}$ space to $f_{j}$ space.

A simple ‘‘projective’’ example with contributing subspace was considered in
Eq. (27) above. The matrix $K^{(\mathbf{f\to x})}_{kk^{\prime}}$ from (25) has
the rank at most $m$ and the best what can be obtained in the projective
paradigm is a solution[13] of ‘‘direct projection’’ type where the least
squares expansion $\mathbf{f}_{LS}(\mathbf{x})$ of $\Ket{f_{j}}$ in
$\Ket{x_{k}}$ (7) is used as the localization point in (13) to obtain the
state $\Ket{\psi_{\mathbf{f}_{LS}(\mathbf{x})}}$ to be used in calculation of
probabilities.

Properly normalized objective function (64) maximizes the probability
transferred from $\mathbf{x}$ to $\mathbf{f}$. Consider a much simpler
problem: find a subspace of $\mathbf{x}$ contributing to the coverage of
$\mathbf{f}$. The $\mathbf{f}$–coverage is determined by
$\mathbf{f}$–Christoffel function $K^{(\mathbf{f})}(\mathbf{g})$ from (23).
Consider it’s values in a $\psi(\mathbf{x})$ state

$\displaystyle\mathrm{Coverage}_{\psi}$
$\displaystyle=\frac{\Braket{\psi}{K^{(\mathbf{f})}}{\psi}}{\Braket{\psi}{\psi}}$
(129)

Previously we considered a similar problem where the Christoffel function $K$
and $\psi$ both were functions on $\mathbf{x}$, see [18], _Appendix B: On The
Christoffel Function Spectrum_. Now the Christoffel function is a function on
$\mathbf{f}$, and $\psi$ is a function on $\mathbf{x}$. The (129) can be
similarly expanded in spectrum of $\mathbf{f}$–Christoffel function matrix

$\displaystyle\Braket{x_{k}}{K^{(\mathbf{f})}}{x_{k^{\prime}}}=\sum\limits_{l=1}^{M}\frac{x^{(l)}_{k}x^{(l)}_{k^{\prime}}}{\sum\limits_{j,j^{\prime}=0}^{m-1}f^{(l)}_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}f^{(l)}_{j^{\prime}}}\omega^{(l)}$
(130)

It is different from (24) with $x$–moments instead of $f$–moments. Consider
generalized eigenvalue problem

$\displaystyle\phi^{[i]}$
$\displaystyle=\sum_{k=0}^{n-1}\alpha^{\phi;[i]}_{k}x_{k}$ $\displaystyle
i=0\dots n-1;$ (131) $\displaystyle\delta_{ii^{\prime}}$
$\displaystyle=\Braket{\phi^{[i]}}{\phi^{[i^{\prime}]}}=\sum_{k,k^{\prime}=0}^{n-1}\alpha^{\phi;[i]}_{k}\Braket{x_{k}x_{k^{\prime}}}\alpha^{\phi;[i^{\prime}]}_{k^{\prime}}$
(132) $\displaystyle\lambda^{[i]}\delta_{ii^{\prime}}$
$\displaystyle=\Braket{\phi^{[i]}}{K^{(\mathbf{f})}}{\phi^{[i^{\prime}]}}=\sum_{k,k^{\prime}=0}^{n-1}\alpha^{\phi;[i]}_{k}\Braket{x_{k}}{K^{(\mathbf{f})}}{x_{k^{\prime}}}\alpha^{\phi;[i^{\prime}]}_{k^{\prime}}$
(133) $\displaystyle\sum\limits_{k^{\prime}=0}^{n-1}$
$\displaystyle\Braket{x_{k}}{K^{(\mathbf{f})}}{x_{k^{\prime}}}\alpha^{\phi;[i]}_{k^{\prime}}=\lambda^{[i]}\sum\limits_{k^{\prime}=0}^{n-1}\Braket{x_{k}x_{k^{\prime}}}\alpha^{\phi;[i]}_{k^{\prime}}$
(134)

Because $\mathbf{x}$– and $\mathbf{f}$– bases are different the condition[18]
$\Braket{1}=\sum_{i=0}^{n-1}\lambda^{[i]}$ no longer holds, it is typically
$\Braket{1}\leq\sum_{i=0}^{n-1}\lambda^{[i]}$ since $m\leq n$; moreover the
sum of $m$ maximal eigenvalues can possibly exceed the total weight
$\Braket{1}\lesseqgtr\sum_{i=0}^{m-1}\lambda^{[i]}$. From Christoffel function
invariance it immediately follows that the sum of $m$ maximal eigenvalues is
equal to $\Braket{1}$ if $\mathbf{f}$ and $\mathbf{x}$ belong to the same
space.

The $m$ eigenstates of (134) corresponding to $m$ maximal eigenvalues
$\lambda^{[i]}$, $i=0\dots n-1$ form the $m$ states contributing most to the
coverage. This is an alternative option for the contributing subspace. The
problem is now reduced to finding a unitary (not partially unitary) operator
$\mathcal{U}$ of the dimension $m\times m$ mapping from $\phi_{k}$ to $f_{j}$,
where $\Ket{\phi}=\sum_{k=0}^{m-1}\phi_{k}\Ket{\phi^{[k]}}$,

$\displaystyle f_{j}$ $\displaystyle=\sum\limits_{k=0}^{m-1}u_{jk}\phi_{k}$
(135)

In this form the optimization problem is greatly simplified and the
$\mathbf{x}$–normalizing in (122) becomes $u_{jk}$ independent:

$\displaystyle\Braket{\psi_{\mathbf{g}}}{\mathcal{U}}{\psi_{\mathbf{\phi}}}^{2}$
$\displaystyle=\frac{\left|\sum\limits_{k=0}^{m-1}\sum\limits_{j,s=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{js}u_{sk}\phi_{k}\right|^{2}}{\sum\limits_{j,j^{\prime}=0}^{m-1}g_{j}G^{\mathbf{f};\,-1}_{jj^{\prime}}g_{j^{\prime}}\sum\limits_{j,j^{\prime}=0}^{m-1}\phi_{j}G^{\phi;\,-1}_{jj^{\prime}}\phi_{j^{\prime}}}$
(136)

This probability is exactly the same as the one we considered above, but with
the $\phi_{k}$ used as the input instead of the $x_{k}$; we also have $n=m$
thus the operator $u_{jk}$ is unitary!

##  Appendix D Software description

  * •

Install java 19 or later.

  * •

Download the latest version of the source code
code_polynomials_quadratures.zip from [20] or from alternative location.

  * •

Decompress and recompile the program. Run a selftest.

        unzip code_polynomials_quadratures.zip
    javac -g com/polytechnik/*/*java
    java com/polytechnik/utils/TestKGO

  * •

Run the program with bundled deterministic data file, test trivial mapping.

        java com/polytechnik/utils/KGO --data_cols=9:0,6:0,4:8:1 \
      --SKtype=FXFX_F_CHRISTOFFEL \
      --approximation=MAXEV_EVADJ \
      --data_file_to_build_model_from=dataexamples/runge_function.csv \
      --output_files_prefix=/tmp/out_

  * •

There are a number of `--approximation=` available options. There are no
perfect implementation yet available.

## References

  * Malyshkin  _et al._ [2001] V. G. Malyshkin, R. Bakhramov, and A. E. Gorodetsky, A Massive Local Rules Search Approach to the Classification Problem, arXiv:cs/0609007 10.48550/arXiv.cs/0609007 (2001), cs/0609007 .
  * Rosenblatt [1958] F. Rosenblatt, The perceptron: a probabilistic model for information storage and organization in the brain., Psychological review 65, 386 (1958).
  * Vapnik and Chervonenkis [1974a] V. Vapnik and A. Y. Chervonenkis, The method of ordered risk minimization, I, Avtomatika i Telemekhanika 8, 21 (1974a).
  * Vapnik and Chervonenkis [1974b] V. Vapnik and A. Y. Chervonenkis, On the method of ordered risk minimization, II, Avtomatika i Telemekhanika 9, 29 (1974b).
  * Hájek and Havránek [1977] P. Hájek and T. Havránek, On generation of inductive hypotheses, International Journal of Man-Machine Studies 9, 415 (1977).
  * Vapnik [2013] V. Vapnik,  _The nature of statistical learning theory_ (Springer science & business media, 2013).
  * Witten and Frank [2002] I. H. Witten and E. Frank, Data mining: practical machine learning tools and techniques with Java implementations, Acm Sigmod Record 31, 76 (2002).
  * Zadeh [1965] L. A. Zadeh, Fuzzy sets, Information and control 8, 338 (1965).
  * Hájek [1995] P. Hájek, Fuzzy logic and arithmetical hierarchy, Fuzzy sets and Systems 73, 359 (1995).
  * Bengio  _et al._ [2013] Y. Bengio, A. Courville, and P. Vincent, Representation learning: A review and new perspectives, IEEE transactions on pattern analysis and machine intelligence 35, 1798 (2013).
  * Malyshkin [2015] V. G. Malyshkin, Norm-Free Radon-Nikodym Approach to Machine Learning, ArXiv e-prints 10.48550/arXiv.1512.03219 (2015), http://arxiv.org/abs/1512.03219, arXiv:1512.03219 [cs.LG] .
  * Bobyl  _et al._ [2020] A. Bobyl, V. Davydov, and V. Malyshkin, On The Radon–Nikodym Machine Learning Parallelization, in  _The 4th International Conference on Future Networks and Distributed Systems (ICFNDS)_ (2020) pp. 1–5.
  * Malyshkin [2019a] V. G. Malyshkin, On The Radon–Nikodym Spectral Approach With Optimal Clustering, arXiv preprint arXiv:1906.00460 10.48550/arXiv.1906.00460 (2019a), arXiv:1906.00460 [cs.LG] .
  * Marx  _et al._ [2019] S. Marx, E. Pauwels, T. Weisser, D. Henrion, and J.-B. Lasserre, Tractable semi-algebraic approximation using Christoffel-Darboux kernel, arXiv preprint arXiv:1904.01833 10.48550/arXiv.1904.01833 (2019).
  * Malyshkin [2019b] V. G. Malyshkin, Market Dynamics: On Directional Information Derived From (Time, Execution Price, Shares Traded) Transaction Sequences., ArXiv e-prints 10.48550/arXiv.1903.11530 (2019b), arXiv:1903.11530 [q-fin.TR] .
  * Shor [2013] N. Z. Shor,  _Nondifferentiable optimization and polynomial problems_, Vol. 24 (Springer Science & Business Media, 2013).
  * Malyshkin [2017] G. S. Malyshkin, The comparative efficiency of classical and fast projection algorithms in the resolution of weak hydroacoustic signals (Сравнительная эффективность классических и быстрых проекционных алгоритмов при разрешении слабых гидроакустических сигналов), Acoustical Physics 63, 216 (2017), doi:10.1134/S1063771017020099 (eng) ; doi:10.7868/S0320791917020095 (рус).
  * Malyshkin [2018] V. G. Malyshkin, On Lebesgue Integral Quadrature, ArXiv e-prints 10.48550/arXiv.1807.06007 (2018), arXiv:1807.06007 [math.NA] .
  * Kraus [1983] K. Kraus,  _States, Effects, and Operations: Fundamental Notions of Quantum Theory_, Lecture Notes in Physics, Vol. 190 (Springer-Verlag, 1983) Lectures in Mathematical Physics at the University of Texas at Austin.
  * Malyshkin [2014] V. G. Malyshkin,  (2014), the code for polynomials calculation, http://www.ioffe.ru/LNEPS/malyshkin/code.html, there is also an alternative location.