Membership-Mappings for Practical Secure Distributed Deep Learning
Mohit Kumar,
Weiping Zhang,
Lukas Fischer,
and Bernhard Freudenthaler
This work was partly supported by the Austrian Research Promotion Agency (FFG) Project SMiLe (Secure Machine Learning Applications with Homomorphically Encrypted Data); FFG COMET-Modul S3AI (Security and Safety for Shared Artificial Intelligence); FFG Sub-Project PETAI (Privacy Secured Explainable and Transferable AI for Healthcare Systems); EU Horizon 2020 Project SERUMS (Securing Medical Data in Smart Patient-Centric Healthcare Systems); the BMK; the BMDW; and the Province of Upper Austria in the frame of the COMET Programme managed by FFG.
M. Kumar is with the Software Competence Center Hagenberg, Austria and the Faculty of Computer Science and Electrical Engineering, University of Rostock, Germany.
W. Zhang is with the Northwestern Polytechnical University, Suzhouzhong Road, Suzhou, China, e-mail<EMAIL_ADDRESS>L. Fischer and B. Freudenthaler are with the Software Competence Center Hagenberg, Austria.

This study leverages the data representation capability of fuzzy based membership-mappings for practical secure distributed deep learning using fully homomorphic encryption. The impracticality issue of secure machine (deep) learning with fully homomorphic encrypted data, arising from large computational overhead, is addressed via applying fuzzy attributes. Fuzzy attributes are induced by globally convergent and robust variational membership-mappings based local deep models. Fuzzy attributes combine the local deep models in a robust and flexible manner such that the global model can be evaluated homomorphically in an efficient manner using a boolean circuit composed of bootstrapped binary gates. The proposed method, while preserving privacy in a distributed learning scenario, remains accurate, practical, and scalable. The method is evaluated through numerous experiments including demonstrations through MNIST dataset and Freiburg Groceries Dataset. Further, a biomedical application related to mental stress detection on individuals is considered.

Membership-mappings, fully homomorphic encryption, fuzzy attributes, distributed deep learning, privacy.

§ INTRODUCTION

Fully homomorphic encryption (FHE) is a solution to the privacy concerns in the cloud computing scenario. The first FHE scheme <cit.> is based on ideal lattices and the bootstrapping procedure is introduced to reduce the noise contained in a ciphertext for allowing arbitrary computations. The bootstrapping operation is performed on a ciphertext via evaluating the decryption function homomorphically using the bootstrapping key (which is the encryption of the private decryption key under the public encryption key). Bootstrapping is the computationally most expensive part of a homomorphic encryption scheme. The theoretical breakthrough of <cit.> was followed by several attempts to develop more practical FHE schemes. The scheme introduced in <cit.> uses only elementary modulo arithmetic and is homomorphic with regard to both addition and multiplication. This scheme was improved in <cit.> with reduced public key size, extended in <cit.> to support encrypting and homomorphically processing a vector of plaintexts as a single ciphertext, and generalized to non-binary messages in <cit.>. Schemes based on a different hard problem, referred to as Learning With Errors (LWE) problem <cit.>, were constructed and many current schemes still rely on LWE or its variants. A FHE scheme constructed in <cit.> is based solely on the standard LWE assumption that is known to be at least as hard as solving hard problems in general lattices. In a variant of the LWE problem, called ring learning with errors problem (RLWE) problem, the algebraic structure of the underlying hard problem reduces the key sizes and speeds up the homomorphic operations. A leveled fully homomorphic encryption scheme based on LWE or RLWE, without bootstrapping procedure, was proposed in <cit.>. The ciphertexts contain a certain amount of noise for security purposes that grows with homomorphic operations. For a better management of the noise growth, <cit.> introduced a modulus switching technique where a complete ladder of moduli is used for scaling down the ciphertext to the next modulus after each multiplication. A tensoring technique for LWE-based FHE that reduced ciphertext noise growth after multiplication from quadratic to linear was introduced in <cit.>. As the scheme of <cit.> does no longer require the rescaling of the ciphertext, this scheme was called a scale-invariant fully homomorphic encryption scheme. An RLWE version of the scale-invariant scheme of <cit.> was created in <cit.>. A technique for building LWE based FHE scheme called as approximate eigenvector method in which homomorphic addition and multiplication are just matrix addition and multiplication was proposed in <cit.>. The essence of this scheme is that the secret key is an approximate eigenvector of the ciphertext matrix and the message is the corresponding eigenvalue. Several works that followed the theoretical breakthrough of <cit.> were aimed at improving the bootstrapping as the bootstrapping remained the bottleneck for an efficient FHE in practice. A much faster bootstrapping, based on a scheme similar to the type of <cit.> that allows to homomorphically compute simple bit operations and refresh (bootstrap) the resulting output in less than a second, was devised in <cit.>. Finally, the TFHE scheme was proposed in <cit.> that features an improved bootstrapping procedure that is considerably more efficient than the previous state of the art. The TFHE scheme generalizes previous structures and schemes over the torus (i.e., the reals modulo 1) and improves the bootstrapping dramatically. For practical applications, TFHE is an open-source C/C++ library <cit.> implementing the ring-variant of <cit.> together with the optimizations of <cit.>. TFHE library implements a very fast gate-by-gate bootstrapping and supports the homomorphic evaluation of the binary gates. The library allows to evaluate homomorphically an arbitrary boolean circuit composed of binary gates without restriction on the number of gates or on their composition, over encrypted data, without decrypting. However, the bootstrapped bit operations are still several times slower than their plaintext equivalents.

The fuzzy systems' capability of handling uncertainties in a rigorous mathematical manner has motivated combining fuzzy theory with deep models <cit.>. Deep neural networks outperform classical machine learning techniques in a wide range of applications but their training requires a large amount of data. The issues, such as determining the optimal model structure, requirement of large training dataset, and iterative time-consuming nature of numerical learning algorithms, are inherent to the neural networks based parametric deep models. The nonparametric approach on the other hand can be promising to address the issue of optimal choice of model structure. However, an analytical solution instead of iterative gradient-based numerical algorithms will be still desired for the learning of deep models. These motivations have led to the development of a nonparametric deep model <cit.> that is learned analytically for representing data points. The study in <cit.> introduces the concept of fuzzy-mapping which is about representing mappings through a fuzzy set such that the dimension of membership function increases with an increasing data size. A relevant result is that a deep autoencoder model formed via a composition of finite number of nonparametric fuzzy-mappings can be learned analytically via variational optimization technique. However, <cit.> didn't provide a formal mathematical framework for the conceptualization of so-called fuzzy-mapping. The study in <cit.> provides to fuzzy-mapping a measure-theoretic conceptualization and refers it to as membership-mapping. Further, the membership-mapping could serve as the building block of deep models <cit.>. An alternative idea of deep autoencoder, that consists of layers such that each layer learns data representation at certain abstraction level through a membership-mappings based autoencoder, is introduced in <cit.> for data representation learning.

The aim of this study is to develop a methodology for practical secure distributed deep learning using fully homomorphic encryption. The machine (deep) learning with fully homomorphically encrypted data remains impractical due to the large computational overhead. Thus, we address in this study the impracticality issue of secure distributed deep learning via

* leveraging data representation learning capability of globally convergent and robust membership-mappings to build local deep models,

* using local deep models to induce fuzzy attributes such that defined fuzzy attributes learn data representation,

* combining local deep models in a robust and flexible manner by means of fuzzy attributes and fuzzy rules to define a global model,

* defining the global model in such a way that the global model can be evaluated homomorphically in an efficient manner using a boolean circuit composed of bootstrapped binary gates,

* implementing very fast gate-by-gate bootstrapping to homomorphically evaluate the global model (that combines the distributed local models) to predict the output.

The proposed approach to secure distributed deep learning is novel. To the best knowledge of the authors, this is the first study to apply fuzzy attributes, which are induced by globally convergent and robust variational membership-mappings based local deep models, for an efficient homomorphic evaluation of the global model. The idea of using fuzzy sets and fuzzy rules to aggregate the local private deep models for building the global model was also considered in <cit.>, however, under differential privacy framework. Differential privacy preserves the privacy of the training dataset via adding random noise to ensure that an adversary can not infer any single data instance by observing model parameters or model outputs. The amount of noise depends upon the value of privacy-loss bound. A major limitation of the differential privacy is that a sufficiently low value of privacy-loss bound results in a loss of accuracy. Moreover, it is not clear how to practically choose the value of privacy-loss bound. FHE approach on the other hand does not lead to the loss of accuracy, however, requires a large computational time.

The data representation learning capability of membership-mappings is central to our methodology. Although <cit.> provided an algorithm for the variational learning of membership-mappings via following the approach of <cit.>, there remains the following two limitations:

* there is no mathematical proof regarding the convergence of the learning algorithm, and

* there is no mathematical analysis on the robustness of the learning algorithm.

This study addresses these two limitations and presents a more simple and elegant estimation algorithm for the variational learning of membership-mappings. A convergence analysis is carried out via deriving a sufficient condition for the convergence of the estimation algorithm. The convergence analysis allows to provide a globally convergent algorithm for the variational learning of membership-mappings based deep models. Further, it is shown that the learning algorithm provides a robust estimation of model parameters via solving a min-max estimation problem. The proposed method for secure distributed deep learning is implemented using MATLAB R2017b and TFHE C/C++ library <cit.>. Experiments have been performed to evaluate the method (in-terms of accuracy and computational time) on MNIST dataset, Freiburg Groceries Dataset, and a biomedical dataset consisting of heart rate interval measurements of different subjects. Further, the scalability of the method as the number of parties participating in collaborative learning increases is investigated.

The paper is organized into sections. Section <ref> reviews the membership-mappings from previous works. An estimation algorithm for the variational learning of membership-mappings is provided in Section <ref>. The convergence and robustness issues have been addressed in Section <ref>. Section <ref> considers the application of membership-mappings to the secure distributed deep learning problem. The experimental validation of the method is provided in Section <ref>. Finally, the concluding remarks are stated in Section <ref>.

§ REVIEW OF MEMBERSHIP-MAPPINGS

This section is dedicated to the review of variational membership-mappings from previous works <cit.>.

§.§ Notations and Definitions

* Let $n,N,p,M \in \mathbb{N}$.

* Let $\mathcal{B}(\mathbb{R}^N)$ denote the Borel $\sigma-$algebra on $\mathbb{R}^N$, and let $\lambda^N$ denote the Lebesgue measure on $\mathcal{B}(\mathbb{R}^N)$.

* Let $(\mathcal{X}, \mathcal{A} , \rho)$ be a probability space with unknown probability measure $\rho$.

* Let us denote by $\mathcal{S}$ the set of finite samples of data points drawn i.i.d. from $\rho$, i.e.,
𝒮     :=      { (x^i  ∼ρ)_i=1^N   |   N ∈ℕ }.

* For a sequence $\mathrm{x} = (x^1,\cdots,x^N) \in \mathcal{S}$, let $|\mathrm{x}|$ denote the cardinality i.e. $|\mathrm{x}| = N$.

* If $\mathrm{x} = (x^1, \cdots, x^N),\; \mathrm{a}= (a^1, \cdots, a^M) \in \mathcal{S}$, then $\mathrm{x} \wedge \mathrm{a}$ denotes the concatenation of the sequences $\mathrm{x}$ and $\mathrm{a}$, i.e., $\mathrm{x} \wedge \mathrm{a} = (x^1, \ldots, x^N, a^1, \ldots, a^M)$.

* Let us denote by $\mathbb{F}(\mathcal{X})$ the set of $\mathcal{A}$-$\mathcal{B}(\mathbb{R})$ measurable functions $f:\mathcal{X} \rightarrow \mathbb{R}$, i.e.,
𝔽(𝒳)     :=      { f:𝒳 →ℝ    |    }.

* For convenience, the values of a function $f \in \mathbb{F}(\mathcal{X})$ at points in the collection $\mathrm{x} = (x^1,\cdots,x^N)$ are represented as $f(\mathrm{x})=(f(x^1),\cdots,f(x^N))$.

* Given two $\mathcal{B}(\mathbb{R}^N)-\mathcal{B}(\mathbb{R})$ measurable mappings, $g:\mathbb{R}^N \rightarrow \mathbb{R}$ and $\mu:\mathbb{R}^N \rightarrow \mathbb{R}$, the weighted average of $g(\mathrm{y})$ over all $\mathrm{y} \in \mathbb{R}^{N}$, with $\mu(\mathrm{y})$ as the weighting function, is computed as
< g >_μ    :=     1/ ∫_ℝ^N  μ(y)  λ^N(y) ∫_ℝ^N g(y) μ(y)  λ^N(y).

* For a sequence $\mathrm{x} \in \mathcal{S}$, assume that a membership function $\zeta_{\mathrm{x}}:\mathbb{R}^{|\mathrm{x}|} \rightarrow [0,1]$ satisfies the following properties:

* $\zeta_{\mathrm{x}}(\mathrm{y}) > 0$ for all $\mathrm{y} \in \mathbb{R}^{|\mathrm{x}|}$, i.e.,

[ζ_x]     =     ℝ^|x|.

* the functions $\zeta_{\mathrm{x}}$ are absolutely continuous and Lebesgue integrable over the whole domain such that for all $\mathrm{x}\in \mathcal{S}$ we have

\begin{eqnarray}
\label{eq:positive}
0 < \int_{\mathbb{R}^{|\mathrm{x}|}} \zeta_{\mathrm{x}}\, \dd \lambda^{|\mathrm{x}|}  < \infty.
\end{eqnarray}

* the membership function induced probability measures $\mathbb{P}_{\zeta_{\mathrm{x}}}$, defined on any $A \in \mathcal{B}(\mathbb{R}^{|\mathrm{x}|})$, as
ℙ_ζ_x(A)     :=      1/ ∫_ℝ^|x| ζ_x  λ^|x| ∫_A ζ_x  λ^|x|

are consistent in the sense that for all $\mathrm{x},\;\mathrm{a} \in \mathcal{S}$:
ℙ_ζ_x ∧a(A ×ℝ^|a|)     =     ℙ_ζ_x(A).

The collection of membership functions satisfying aforementioned assumptions is denoted by
Θ    :=     { ζ_x:ℝ^|x| →[0,1]   |   (<ref>),  (<ref>), (<ref>),  x ∈𝒮}.

A Student-t membership-mapping, $\mathcal{F} \in \mathbb{F}(\mathcal{X})$, is a mapping with input space $\mathcal{X} = \mathbb{R}^n$ and a membership function $\zeta_{\mathrm{x}} \in \Theta$ that is Student-t like:

ζ_x(y)     =     (1 + 1/ν- 2  ( y - m_y )^T K^-1_xx ( y- m_y) )^-ν+|x|/2
where $\mathrm{x} \in \mathcal{S}$, $\mathrm{y} \in \mathbb{R}^{|\mathrm{x}|}$, $\nu \in \mathbb{R}_{+}\setminus [0,2]$ is the degrees of freedom, $\mathrm{m}_{\mathrm{y}} \in \mathbb{R}^{|\mathrm{x}|}$ is the mean vector, and $K_{\mathrm{x}\mathrm{x}} \in \mathbb{R}^{|\mathrm{x}| \times |\mathrm{x}|}$ is the covariance matrix with its $(i,j)-$th element given as
(K_xx)_i,j     =     kr(x^i,x^j)

where $kr: \mathbb{R}^n \times \mathbb{R}^n \rightarrow \mathbb{R}$ is a positive definite kernel function defined as
kr(x^i,x^j)     =      σ^2 exp(-0.5∑_k = 1^n w_k |  x^i_k - x^j_k |^2)

where $x_k^i$ is the $k-$th element of $x^i$, $\sigma^2$ is the variance parameter, and $w_{k} \geq 0$ (for $k \in \{1,\cdots,n\}$).

§.§ Conditionally Deep Autoencoders

This subsection reviews the conditionally deep models <cit.> formed by the compositions of membership-mappings.

A membership-mapping autoencoder, $\mathcal{G}:\mathbb{R}^p \rightarrow \mathbb{R}^p$, maps an input vector $y \in \mathbb{R}^p$ to $\mathcal{G}(y) \in \mathbb{R}^p$ such that
𝒢(y)        def=      [[][c],c/c/c,  ℱ_1(Py)     ⋯     ℱ_p(Py)

where $\mathcal{F}_j$ ($j \in \{1,2,\cdots,p\}$) is a Student-t membership-mapping, $P \in \mathbb{R}^{n \times p} (n \leq p)$ is a matrix such that the product $Py$ is a lower-dimensional encoding for $y$.

A conditionally deep membership-mapping autoencoder, $\mathcal{D}:\mathbb{R}^p \rightarrow \mathbb{R}^p$, maps a vector $y \in \mathbb{R}^p$ to $\mathcal{D}(y) \in \mathbb{R}^p$ through a nested composition of finite number of membership-mapping autoencoders such that
y^l     =     (𝒢_l ∘⋯∘𝒢_2∘𝒢_1)(y),   ∀l ∈{1,2,⋯,L }

l^*     =      min_l   ∈   {1,2,⋯,L }   y -  y^l ^2

𝒟(y)     =     y^l^*,

where $\mathcal{G}_l(\cdot)$ is a membership-mapping autoencoder (Definition <ref>).

A wide CDMMA, $\mathcal{WD}:\mathbb{R}^p \rightarrow \mathbb{R}^p$, maps a vector $y \in \mathbb{R}^p$ to $\mathcal{WD}(y) \in \mathbb{R}^p$ through a parallel composition of $S$  ($S \in \mathcal{Z}_+$) number of CDMMAs such that
𝒲𝒟(y)     =     𝒟_s^*(y)

s^*     =      min_s ∈{1,2,⋯,S }   y -  𝒟_s(y)  ^2,

where $\mathcal{D}_s(y)$ is the output of $s-$th CDMMA.

§ VARIATIONAL LEARNING OF MEMBERSHIP-MAPPINGS

§.§ A Modeling Scenario

Given a dataset $\{(x^i,y^i)\;|\; x^i \in \mathbb{R}^n,\;y^i \in \mathbb{R}^p,\; i \in \{1,\cdots,N \} \}$, it is assumed that there exist zero-mean Student-t membership-mappings $\mathcal{F}_1, \cdots, \mathcal{F}_p \in \mathbb{F}(\mathbb{R}^n)$ such that
y^i    ≈    [[][c],c/c/c, ℱ_1(x^i)      ⋯    ℱ_p(x^i)  ]^T.

For $ j \in \{1,2,\cdots,p\}$, define
y_j     =     [[][c],c/c/c,y_j^1     ⋯   y_j^N ]^T ∈ℝ^N

f_j     =     [[][c],c/c/c,ℱ_j(x^1)     ⋯    ℱ_j(x^N) ]^T ∈ℝ^N

where $y_j^i$ denotes the $j-$th element of $y^i$. A set of auxiliary inducing points, $\mathrm{a}  =  \{a^{m} \in \mathbb{R}^{n}\; | \; m \in\{1,\cdots,M\} \}$, is introduced to define
u_j =  [[][c],c/c/c,ℱ_j(a^1)     ⋯    ℱ_j(a^M) ]^T ∈ℝ^M.

§.§ Membership Functional Representation Approach

It follows from <cit.> that $\mathrm{f}_j$, based upon an interpolation on the auxiliary-outputs $\mathrm{u}_j $, is represented by means of a membership function, $\mu_{\mathrm{f}_j;\mathrm{u}_j}:\mathbb{R}^N \rightarrow [0,1]$, as
( μ_f_j;u_j(f̃_j) )^-2/ν+M+N  =  1 +

(f̃_j -  m̅_f_j)^T  ( ν+ (u_j)^T K_aa^-1 u_j  - 2/ν+ M - 2 K̅_xx )^-1(f̃_j -  m̅_f_j)/ν+ M - 2

m̅_f_j     =      K_xa K_aa^-1  u_j

K̅_xx     =     K_xx - K_xa K_aa^-1 K_xa^T,

where $K_{\mathrm{a}\mathrm{a}} \in \mathbb{R}^{M \times M}$ and $K_{\mathrm{x}\mathrm{a}}  \in \mathbb{R}^{N \times M}$ are matrices with their $(i,j)-$th elements given as
( K_aa )_i,j     =     kr(a^i,a^j)

( K_xa )_i,j     =     kr(x^i,a^j)

where $kr: \mathbb{R}^n \times \mathbb{R}^n \rightarrow \mathbb{R}$ is a positive definite kernel function defined as in (<ref>).

$\mathrm{y}_j$, for given $\mathrm{f}_j$, is represented by means of a membership function, $\mu_{\mathrm{y}_j;\mathrm{f}_j}:\mathbb{R}^N \rightarrow [0,1]$, as
μ_y_j;f_j(ỹ_j)     =     exp(- 0.5 β  ỹ_j -  f_j ^2 )

where $\beta > 0$ is the precision value.

$\mathrm{y}_j$, for given $\mathrm{u}_j$, is represented by means of a membership function, $\mu_{\mathrm{y}_j;\mathrm{u}_j}:\mathbb{R}^N \rightarrow [0,1]$, as
μ_y_j;u_j(ỹ_j)      ∝    exp( < log( μ_y_j;f_j(ỹ_j) )  >_μ_f_j;u_j )
where $\mu_{\mathrm{y}_j;\mathrm{f}_j}$ is given by (<ref>), $\mu_{\mathrm{f}_j;\mathrm{u}_j}$ is defined as in (<ref>), and $<\cdot>_{\cdot}$ is the averaging operation as defined in (<ref>). Thus, $\mu_{\mathrm{y}_j;\mathrm{u}_j}$ is obtained from $\log( \mu_{\mathrm{y}_j;\mathrm{f}_j})$ after averaging out the variables  $\mathrm{f}_j$ using its membership function. It can be shown that
μ_y_j;u_j(ỹ_j)  ∝exp(  - 0.5β ỹ_j^2  + (u_j)^T K̂_u_j^-1 m̂_u_j (ỹ_j) .

. - 0.5 (u_j)^T K̂_u_j^-1 u_j + 0.5 (u_j)^T  K_aa^-1  u_j +  {/ (ỹ_j,u_j)} )
where $\hat{K}_{\mathrm{u}_j}$, $\hat{m}_{\mathrm{u}_j}(\tilde{\mathrm{y}}_j)$ are given as
(K̂_u_j)^-1     =     K_aa^-1  + βK_aa^-1  K_xa^T K_xa  K_aa^-1

+  βtr(K_xx - K_aa^-1  K_xa^T K_xa )/ν+ M - 2 K_aa^-1

m̂_u_j(ỹ_j)     =     βK̂_u_j K_aa^-1 K_xa^T ỹ_j,

and $ \{/ (\tilde{\mathrm{y}}_j,\mathrm{u}_j)\}$ represents all those terms which are independent of both $\tilde{\mathrm{y}}_j$ and $\mathrm{u}_j$. The constant of proportionality in (<ref>) is chosen to exclude $ (\tilde{\mathrm{y}}_j,\mathrm{u}_j)-$independent terms in the expression for $\mu_{\mathrm{y}_j;\mathrm{u}_j}$, i.e.,
μ_y_j;u_j(ỹ_j)   =  exp(  (u_j)^T K̂_u_j^-1 m̂_u_j (ỹ_j) .

.  - 0.5 (u_j)^T K̂_u_j^-1 u_j + 0.5 (u_j)^T  K_aa^-1  u_j - 0.5β ỹ_j^2 ).

$\mathrm{y}_j$ is represented by means of a membership function, $\mu_{\mathrm{y}_j} : \mathbb{R}^N \rightarrow [0,1]$, as
μ_y_j(ỹ_j)     ∝    exp( <  log( μ_y_j;u_j (ỹ_j) )>_μ_u_j )

where $ \mu_{\mathrm{y}_j;\mathrm{u}_j}$ is given by (<ref>) and $\mu_{\mathrm{u}_j}:  \mathbb{R}^M \rightarrow [0,1]$ is a membership function representing $\mathrm{u}_j$. Thus, $\mu_{\mathrm{y}_j}$ is obtained from $\log(\mu_{\mathrm{y}_j;\mathrm{u}_j})$ after averaging out the auxiliary-outputs $\mathrm{u}_j$ using membership function $\mu_{\mathrm{u}_j}$.

§.§ Variational Optimization of Data-Model

The data model (<ref>) involves the membership function $\mu_{\mathrm{u}_j}$. To determine $\mu_{\mathrm{u}_j}$ for a given $\mathrm{y}_j$, $\log(\mu_{\mathrm{y}_j}(\mathrm{y}_j))$ is maximized w.r.t. $\mu_{\mathrm{u}_j}$ around an initial guess. The zero-mean Gaussian membership function with covariance as equal to $K_{\mathrm{a}\mathrm{a}}$ is taken as the initial guess. It follows from (<ref>) that maximization of $\log(\mu_{\mathrm{y}_j}(\mathrm{y}_j))$ is equivalent to the maximization of $\left<  \log( \mu_{\mathrm{y}_j;\mathrm{u}_j}(\mathrm{y}_j) )\right>_{\mu_{\mathrm{u}_j}}$.

The solution of following maximization problem:
μ^*_u_j       =        max_μ_u_j   [ <  log( μ_y_j;u_j(y_j) )>_μ_u_j .

. - < log( μ_u_j(u_j)/exp(-0.5 (u_j)^T K_aa^-1 u_j  ) ) >_ μ_u_j ]

under the fixed integral constraint:
∫_ℝ^M μ_u_j    λ^M   =  C_u_j > 0

where the value of $C_{\mathrm{u}_j}$ is so chosen such that the maximum possible values of $\mu^*_{\mathrm{u}_j} $ remain as equal to unity, is given as
μ^*_u_j(u_j)   =

exp(- 0.5 (u_j - m̂_u_j(y_j))^T K̂_u_j^-1(u_j - m̂_u_j(y_j)) )

where $\hat{K}_{\mathrm{u}_j}$ and $\hat{m}_{\mathrm{u}_j}$ are given by (<ref>) and (<ref>) respectively. The solution of the optimization problem results in
μ_y_j(ỹ_j)   ∝exp({/ (y_j,ỹ_j)}   .

. -  0.5β{ ỹ_j^2  - 2 (m̂_u_j (y_j))^T K_aa^-1 K_xa^T ỹ_j . .

. .  +   (m̂_u_j (y_j))^TK_aa^-1  K_xa^T K_xa K_aa^-1m̂_u_j (y_j)  . .

.  .  +  tr(K_xx - K_aa^-1  K_xa^T K_xa )/ν+ M - 2  (m̂_u_j (y_j))^T K_aa^-1m̂_u_j (y_j)  } )

where $\{/ (\mathrm{y}_j,\tilde{\mathrm{y}}_j)\}$ represents all $(\mathrm{y}_j,\tilde{\mathrm{y}}_j)-$independent terms.

The proof is similar to as that of Result 2 in <cit.>.

The constant of proportionality in (<ref>) is chosen to exclude $(\mathrm{y}_j,\tilde{\mathrm{y}}_j)-$independent terms resulting in
μ_y_j(ỹ_j)  =   exp( - β/2 { ỹ_j^2 - 2 (m̂_u_j (y_j))^T K_aa^-1 K_xa^T ỹ_j . .

. .  +  (m̂_u_j (y_j))^TK_aa^-1  K_xa^T K_xa K_aa^-1m̂_u_j (y_j)  . .

.  .  + (m̂_u_j (y_j))^Ttr(K_xx - K_aa^-1  K_xa^T K_xa )/ν+ M - 2 K_aa^-1m̂_u_j (y_j)  } ).

§.§ Estimation of Membership-Mapping Parameters

$\mathcal{F}_j(x^{i})$ (which is the $i-$th element of vector $\mathrm{f}_j$ (<ref>)) can be estimated as
ℱ_j(x^i)     :=      <  < (f_j)_i >_μ_f_j;u_j    >_ μ^*_u_j

where $(\mathrm{f}_j)_i$ denotes the $i-$th element of $\mathrm{f}_j$, $\mu_{\mathrm{f}_j;\mathrm{u}_j } $ is defined as in (<ref>), and $\mu^*_{\mathrm{u}_j}$ is the optimal membership function (<ref>) representing $\mathrm{u}_j$. That is, $\mathcal{F}_j(x^{i})$, being a function of $\mathrm{u}_j$, is averaged over $\mathrm{u}_j$ for an estimation.

Let $G(x) \in \mathbb{R}^{1 \times M}$ be a vector-valued function defined as
G(x)    :=      [[][c],c/c/c,kr(x,a^1)     ⋯    kr(x,a^M)  ]

where $kr: \mathbb{R}^n \times \mathbb{R}^n \rightarrow \mathbb{R}$ is defined as in (<ref>). It is shown in Appendix <ref> that
ℱ_j(x^i)      =     G(x^i) ( K_xa^T K_xa   +  τK_aa + β^-1 K_aa)^-1  K_xa^T y_j
where $\tau$ is given as
τ    :=     tr(K_xx - K_aa^-1  K_xa^T K_xa )/ν+ M - 2.

Define a vector $\alpha_j \in \mathbb{R}^M$ as
α_j(β^-1)     :=     ( K_xa^T K_xa   + τK_aa   +    β^-1 K_aa )^-1  (K_xa )^T y_j
so that $\widehat{ \mathcal{F}_j(x^{i})}$ could be expressed as
ℱ_j(x^i)     =     (G(x^i) ) α_j(β^-1).

It follows from (<ref>) that estimation of the membership-mapping outputs requires computing $\alpha_j$ via (<ref>) which in-turn requires estimating the inverse precision value $\beta^{-1}$. The inverse precision value $\beta^{-1}$ is iteratively estimated as the inverse of the mean squared error between data and membership-mappings outputs. That is,
β^-1     =     1/pN∑_j=1^p ∑_i=1^N |y_j^i -   ℱ_j(x^i) |^2

where $\widehat{ \mathcal{F}_j(x^{i})}$ is the estimated membership-mapping output given as in (<ref>). Expression (<ref>) using (<ref>) can be expressed as
β^-1     =     1/pN∑_j=1^p ∑_i=1^N |y_j^i -  (G(x^i) ) α_j(β^-1) |^2.

As $G(x^i)$ is equal to the $i-$th row of matrix $K_{\mathrm{x}\mathrm{a}}$, (<ref>) can be expressed as
β^-1     =     1/pN∑_j=1^p  y_j  - K_xa α_j(β^-1) ^2.

We suggest to estimate $\beta^{-1}$ and $\alpha_j$ iteratively using (<ref>) and (<ref>) till the convergence.

§ A GLOBALLY CONVERGENT LEARNING ALGORITHM AND ROBUSTNESS ANALYSIS FOR VARIATIONAL MEMBERSHIP-MAPPINGS

§.§ Convergence Analysis

In this subsection, we study the convergence of estimation algorithm (<ref>, <ref>). In particular, we derive a sufficient condition for estimation algorithm (<ref>, <ref>) to converge. For this, consider the singular value decomposition of $K_{\mathrm{x}\mathrm{a}} $:
K_xa      =     U [  S
0 ] V^T

where $U \in \mathbb{R}^{N \times N}$ and $V \in \mathbb{R}^{M \times M}$ are orthogonal, and $S = \text{diag}(s_1,\cdots,s_M)$ is a diagonal matrix with $s_1 \geq s_2 \geq \cdots \geq s_M \geq 0$ being the singular values of $K_{\mathrm{x}\mathrm{a}} $. The vectors $b_j^1 \in \mathbb{R}^M$ and $b_j^2 \in \mathbb{R}^{N-M}$ are defined as
[  b_j^1
b_j^2 ]     =     U^T y_j .

The expression (<ref>) for $\alpha_j $ can be rewritten as
α_j(β^-1)     =     ( V S^2 V^T   + τK_aa   +  β^-1  K_aa )^-1  V S b_j^1 .

y_j  - K_xa α_j(β^-1)  =   U U^T y_j

-  U [  S
0 ] V^T ( V S^2 V^T   + (τ+  β^-1 ) K_aa )^-1  V S b_j^1.
Using matrix inversion lemma,
y_j  - K_xa α_j(β^-1)=

U [ ( I + 1/τ+ β^-1 S V^T  K_aa^-1 VS   )^-1  b_j^1

b_j^2  ],
and hence,
y_j  - K_xa α_j(β^-1)^2  =    b_j^2 ^2

+  (τ+ β^-1)^2(b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2  b_j^1

Using (<ref>) in (<ref>),
β^-1   =  1/pN∑_j=1^p    b_j^2 ^2

+   (τ+β^-1)^2/pN ∑_j=1^p (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2  b_j^1.

[A Function in $\beta^{-1}$]
Let $\mathcal{R}$ be a function in $\beta^{-1}$ defined as
ℛ(β^-1)  : =  1/pN∑_j=1^p    b_j^2 ^2

+ (τ+ β^-1)^2/pN∑_j=1^p  (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2  b_j^1.

We have followings:

* $\mathcal{R}(\beta^{-1}) \in \left(\beta^{-1}|_{low}, \beta^{-1}|_{up}  \right), \; \forall \beta^{-1} \in \mathbb{R}_{>0}$, where
β^-1|_low      =     1/pN∑_j=1^p    b_j^2 ^2

β^-1|_up      =      1/pN∑_j=1^p   y_j  ^2.

* The lower and upper bounds on the derivative of $\mathcal{R}$ w.r.t. $\beta^{-1}$ are given as

\begin{align}
0  <  \frac{\dd \mathcal{R}(\beta^{-1})}{\dd \beta^{-1}}  <  \frac{2}{pN} \frac{1}{\tau +\beta^{-1}} \sum_{j=1}^p  \|  b_j^1 \|^2. \label{eq_738573.4562}
\end{align}

* $\mathcal{R}(\beta^{-1})$ has at least one fixed point in $\left(\beta^{-1}|_{low}, \beta^{-1}|_{up}  \right)$.

The proof is provided in Appendix <ref>.

If $\tau$ is chosen such that
τ    >     2/pN∑_j=1^p   y_j  ^2 - 1/pN∑_j=1^p    b_j^2 ^2,

then the iterations
β^-1|_it+1     =     ℛ(β^-1|_it), it ∈{0,1,2,⋯},

with $\beta^{-1}|_{0} \in \left(\beta^{-1}|_{low}, \beta^{-1}|_{up}  \right)$, converge to the unique fixed point of $\mathcal{R}(\beta^{-1})$.

The proof is provided in Appendix <ref>.

§.§ Learning Algorithm

Result <ref> allows to design a globally convergent algorithm for the variational learning of membership-mappings via ensuring the sufficient condition (<ref>). For this, it is observed that $\tau$, defined as in (<ref>), is a function of parameters set $\{M, \sigma^2,\{w_i\}_{i=1}^n,\mathrm{x},\mathrm{a},\nu\}$. It follows from the kernel function definition (<ref>) that
τ(M, σ^2,{w_i}_i=1^n,x,a,ν)     =     σ^2τ(M, 1,{w_i}_i=1^n,x,a,ν).
Using (<ref>), the condition (<ref>) can be rewritten as
σ^2     >     1/τ(M, 1,{w_i}_i=1^n,x,a,ν) ∑_j=1^p(2 y_j  ^2-   b_j^2 ^2)/pN.
To hold the inequality (<ref>), the value of $\sigma^2$ is adjusted as in the following:

$\tau(M, 1,\{w_i\}_{i=1}^n,\mathrm{x},\mathrm{a},\nu) > \frac{\sum_{j=1}^p(2\| \mathrm{y}_j  \|^2- \|  b_j^2 \|^2)}{pN}$
$\sigma^2 = 1$
$\sigma^2 =  \frac{1.1 }{\tau(M, 1,\{w_i\}_{i=1}^n,\mathrm{x},\mathrm{a},\nu)} \frac{\sum_{j=1}^p(2\| \mathrm{y}_j  \|^2- \|  b_j^2 \|^2)}{pN} $.

Further, some of the parameters must be prior chosen which are suggested to be chosen as in the following:

Auxiliary inducing points
The auxiliary inducing points are suggested to be chosen as the cluster centroids:
a = { a^m}_m=1^M  =  cluster_centroid(  {x^i}_i=1^N, M)

where $cluster\_centroid(  \{ x^i \}_{i=1}^N,M)$ represents the k-means clustering on $ \{ x^i \}_{i=1}^N$.

Degrees of freedom
The degrees of freedom associated to the Student-t membership-mapping $\nu \in \mathbb{R}_{+} \setminus [0,2]$ is chosen as
ν    =     2.1

Parameters $(w_1,\cdots,w_n)$
The parameters $(w_1,\cdots,w_n)$ for kernel function (<ref>) are chosen such that $w_{k}$ (for $k\in \{1,2,\cdots,n\}$) is given as
w_k     =     (max_1 ≤i ≤N(x^i_k) - min_1 ≤i ≤N(x^i_k))^-2

where $x^i_k$ is the $k-$th element of vector $x^i \in \mathbb{R}^n$.

Finally, Algorithm <ref> presents a systematic procedure for the variational learning of membership-mappings while ensuring the sufficient condition (<ref>) for the convergence.

A globally convergent algorithm for the variational learning of membership-mappings

[1]
Dataset $\left\{ (x^i,y^i) \; | \;  i \in \{1,\cdots,N \} \right \}$ and maximum possible number of auxiliary points $M_{max} \in \mathbb{Z}_+$ with $M_{max} \leq N$.
Choose $\nu$ and $w = (w_1,\cdots,w_n)$ as in (<ref>) and (<ref>) respectively.
Set iteration count $it = 0$, $M|_0 = M_{max}$, and determine $\mathrm{a}|_0 = \{ a^{m}|_0\}_{m=1}^{M|_0}$ using (<ref>).
$\tau(M|_{it},1,\{w_i\}_{i=1}^n,\mathrm{x},\mathrm{a}|_{it},\nu) \leq 0$
$M|_{it+1} =  M|_{it} - 1$
Determine $\mathrm{a}|_{it+1} = \{ a^{m}|_{it+1}\}_{m=1}^{M|_{it+1}}$ using (<ref>).
$it \leftarrow it + 1$
Set $M = M|_{it}$ and compute $\mathrm{a} = \{ a^{m}\}_{m=1}^M$ using (<ref>).
Compute $K_{\mathrm{x}\mathrm{a}}$ using (<ref>) taking $\sigma^2 = 1$ and perform singular value decomposition of $K_{\mathrm{x}\mathrm{a}}$ to compute orthogonal matrix $U$ such that (<ref>) holds. Further compute $b_j^1$ and $b_j^2$ using (<ref>) for all $j \in \{1,\cdots,p\}$.
$\tau(M, 1,\{w_i\}_{i=1}^n,\mathrm{x},\mathrm{a},\nu) > \frac{\sum_{j=1}^p(2\| \mathrm{y}_j  \|^2- \|  b_j^2 \|^2)}{pN}$
$\sigma^2 = 1$
$\sigma^2 =  \frac{1.1 }{\tau(M, 1,\{w_i\}_{i=1}^n,\mathrm{x},\mathrm{a},\nu)} \frac{\sum_{j=1}^p(2\| \mathrm{y}_j  \|^2- \|  b_j^2 \|^2)}{pN} $.
Compute $\mathrm{a} = \{ a^{m}\}_{m=1}^M$ using (<ref>), $K_{\mathrm{x}\mathrm{x}}$ using (<ref>), $K_{\mathrm{a}\mathrm{a}}$ using (<ref>), and $K_{\mathrm{x}\mathrm{a}}$ using (<ref>).
Compute $\tau$ using (<ref>).
Set iteration count $it = 0$ and $\beta^{-1}|_0 = 0.5(\beta^{-1}|_{low} + \beta^{-1}|_{up})$, where $\beta^{-1}|_{low}$ and $\beta^{-1}|_{up}$ are given by (<ref>) and (<ref>) respectively.
Determine the unique fixed point of $\mathcal{R}(\beta^{-1})$, say $\hat{\beta}^{-1}$, using iterations (<ref>).
Compute matrix $\alpha = \left[\begin{IEEEeqnarraybox*}[][c]{,c/c/c,}  \alpha_1(\hat{\beta}^{-1}) & \cdots & \alpha_p(\hat{\beta}^{-1})
\end{IEEEeqnarraybox*} \right] \in \mathbb{R}^{M \times p}$, where $\alpha_j(\hat{\beta}^{-1})$, $j\in \{1,\cdots,p\}$, is computed using (<ref>).
the parameters set $\mathbb{M} = \{\alpha, \mathrm{a}, M,\sigma,w,\hat{\beta}\}$.

Given the parameters set $\mathbb{M} = \{\alpha, \mathrm{a}, M,\sigma,w,\hat{\beta}\}$ returned by Algorithm <ref>, the learned membership-mappings could be used to predict output corresponding to any arbitrary input data point $x \in \mathbb{R}^n$ as
ŷ(x;𝕄)     =     [[][c],c/c/c,  ℱ_1(x)     ⋯     ℱ_p(x) ]^T

where $\widehat{ \mathcal{F}_j(x)}$, defined as in (<ref>), is the estimated output of $j-$th membership-mapping. It follows from (<ref>) that
ŷ(x;𝕄)     =     (α(β̂))^T(G(x))^T

where $G(\cdot) \in \mathbb{R}^{1 \times M}$ is a vector-valued function (<ref>).

As a CDMMA consists of membership-mapping compositions, Algorithm <ref> can be directly applied for their learning as in Algorithm <ref> and Algorithm <ref>. For practical applications, Algorithm <ref> is suggested for the learning of wide CDMMA where a computational optimization of a free parameter is performed via minimizing the estimated variance of the mean squared error between data and membership-mappings outputs.

A globally convergent algorithm for the variational learning of CDMMA

[1]
Data set $\mathbf{Y} = \left\{ y^i \in \mathbb{R}^p \; | \; i \in \{1,\cdots,N \} \right \}$; the subspace dimension $n \in \{1,2,\cdots,p \}$; maximum number of auxiliary points $M_{max} \in \mathbb{Z}_+$ with $M_{max} \leq N$; the number of layers $L \in \mathbb{Z}_{+}$.
$l=1$ to $L$
Set subspace dimension associated to $l-$th layer as $n_l = \max(n - l + 1,1)$.
Define $P^l \in \mathbb{R}^{n_l \times p}$ such that $i-$th row of $P^l$ is equal to transpose of eigenvector corresponding to $i-$th largest eigenvalue of sample covariance matrix of data set $\mathbf{Y} $.
Define a latent variable $x^{l,i} \in \mathbb{R}^{n_l}$, for $i \in \{1,\cdots,N \}$, as
x^l,i    :=    {
P^ly^i     if $l=1$,

P^l ŷ^l-1(x^l-1,i;𝕄^l-1)     if $l > 1$

where $\hat{y}^{l-1}$ is the estimated output of the $(l-1)-$th layer computed using (<ref>) for the parameters set $\mathbb{M}^{l-1} = \{\alpha^{l-1}, \mathrm{a}^{l-1}, M^{l-1}, \sigma^{l-1},w^{l-1}\}$.
Define $M_{max}^l$ as
M_max^l    :=    {
M_max     if $l=1$,

M^l-1     if $l > 1$

Compute parameters set $\mathbb{M}^l = \{\alpha^{l}, \mathrm{a}^{l}, M^{l}, \sigma^{l},w^{l}, \hat{\beta}^{l}\}$, characterizing the membership-mappings associated to $l-$th layer, using Algorithm <ref> on data set $\left\{ (x^{l,i},y^i) \; | \;  i \in \{1,\cdots,N \} \right \}$ with maximum possible number of auxiliary points $M_{max}^l$.
the parameters set $\mathcal{M} = \{\{\mathbb{M}^1,\cdots,\mathbb{M}^L\}, \{P^1,\cdots,P^L \} \}$.

Given a CDMMA with its parameters being represented by a set $\mathcal{M} = \{\{\mathbb{M}^1,\cdots,\mathbb{M}^L\}, \{P^1,\cdots,P^L \} \}$, the autoencoder can be applied for filtering a given input vector $y \in \mathbb{R}^p$ as follows:
x^l(y;ℳ)    =    {
P^ly,     $l=1$

P^l  ŷ^l-1(x^l-1;𝕄^l-1)     $l \geq 2$

Here, $\hat{y}^{l-1}$ is the output of the $(l-1)-$th layer estimated using (<ref>). Finally, CDMMA's output, $\mathcal{D}(y;\mathcal{M})$, is given as
𝒟(y;ℳ)     =      ŷ^l^*(x^l^*;𝕄^l^*)

l^*      =      min_l   ∈   {1,⋯,L }  y - ŷ^l(x^l;𝕄^l) ^2.

A globally convergent algorithm for the variational learning of wide CDMMA

[1]
Data set $\mathbf{Y} = \left\{ y^i \in \mathbb{R}^p \; | \; i \in \{1,\cdots,N \} \right \}$; the number of layers $L \in \mathbb{Z}_{+}$; the subspace dimension $n \in \{1,2,\cdots,p\}$; an array of possible $r_{max}$ values $\{r_{max}^1,\cdots,r_{max}^{N_r}\}$ with $0< r_{max}^1 < r_{max}^2 < \cdots < r_{max}^{N_r} \leq 1$.
Apply k-means clustering to partition $\mathbf{Y} $ into $S$ subsets, $\{\mathbf{Y}^1, \cdots, \mathbf{Y}^S  \}$, where $S = \lceil N/1000 \rceil$.
$s = 1$ to $S$
$r = r_{max}^1$ to $r = r_{max}^{N_r}$
Apply Algorithm <ref> on $\mathbf{Y}^s$ to build a single-layered CDMMA, $\mathcal{M}^{r} = \{\{\mathbb{M}^{1,r}\}, \{P^{1,r} \} \}$ where $\mathbb{M}^{1,r} = \{\alpha^{1,r}, \mathrm{a}^{1,r}, M^{1,r}, \sigma^{1,r},w^{1,r}, \hat{\beta}^{1,r}\}$, taking $n$ as the subspace dimension; maximum number of auxiliary points as equal to $r \times \#\mathbf{Y}^s$ (where $\#\mathbf{Y}^s$ is the number of data points in $\mathbf{Y}^s$); and $L = 1$.
Set $r_{max} = \arg \; \max_{r \in \{r_{max}^1, r_{max}^2, \cdots, r_{max}^{N_r}\}} \hat{\beta}^{1,r}$.
Build a CDMMA, $\mathcal{M}^s$, by applying Algorithm <ref> on $\mathbf{Y}^s$ taking $n$ as the subspace dimension; maximum number of auxiliary points as equal to $r_{max} \times \#\mathbf{Y}^s$ (where $\#\mathbf{Y}^s$ is the number of data points in $\mathbf{Y}^s$); and $L$ as the number of layers.
the parameters set $\mathcal{P} = \{\mathcal{M}^s\}_{s=1}^S$.

Given a wide CDMMA with its parameters being represented by a set $\mathcal{P} = \{\mathcal{M}^s\}_{s=1}^S$, the autoencoder can be applied for filtering a given input vector $y \in \mathbb{R}^p$ as follows:
𝒲𝒟(y;𝒫)     =       𝒟(y;ℳ^s^*)

s^*     =      min_s ∈{1,2,⋯,S }   y -  𝒟(y;ℳ^s)  ^2,

where $ \widehat{\mathcal{D}}(y;\mathcal{M}^{s})$ is the output of $s-$th CDMMA estimated using (<ref>).

§.§ Robustness Analysis

Consider that the data samples $\{ (x^i,y^i) \; | \; i \in \{1,\cdots,N\}\}$ are subject to deterministic perturbations and thus matrix $K_{\mathrm{x}\mathrm{a}}$ (<ref>) and vector $\mathrm{y}_j$ (<ref>) are subject to perturbations. Let $\Delta K_{\mathrm{x}\mathrm{a}} \in \mathbb{R}^{N \times M}$ and $\Delta \mathrm{y}_j \in \mathbb{R}^{N}$ be the unknown (but bounded) perturbations in $K_{\mathrm{x}\mathrm{a}}$ and $\mathrm{y}_j$ respectively due to the perturbations in the data samples. The data model assumes that there exists some parameters vector $\alpha^*_j \in \mathbb{R}^M$ such that
y_j + Δy_j     =     (K_xa + ΔK_xa)α_j^*.

The modeling problem is concerned with the estimation of $\alpha_j^*$ in the presence of unknown perturbations $\Delta K_{\mathrm{x}\mathrm{a}}$ and $\Delta \mathrm{y}_j $. To show that Algorithm <ref> provides a robust estimation of $\alpha_j^*$, define
Δ_x     :=     ΔK_xa (K_aa^1/2)^-1

where $K_{\mathrm{a}\mathrm{a}}^{1/2}$ is the unique square root of positive definite matrix $K_{\mathrm{a}\mathrm{a}} $ (<ref>). The set of equations (<ref>) can be expressed as
y_j + Δy_j     =     (K_xa + Δ_x K_aa^1/2)α_j^*.

The perturbation matrix $\left[\begin{IEEEeqnarraybox*}[][c]{,c/c,} \Delta_x &  \Delta \mathrm{y}_j\end{IEEEeqnarraybox*} \right]$ is unknown, however, is assumed bounded. That is, there exists a scalar $\delta_{m} > 0$ such that $\| \left[\begin{IEEEeqnarraybox*}[][c]{,c/c,} \Delta_x &  \Delta \mathrm{y}_j\end{IEEEeqnarraybox*} \right] \|_F \leq \delta_{m} $, where $\| \cdot \|_F$ denotes the Frobenius norm. A robust solution to the estimation of parameters seeks to alleviate the worst-case effect of perturbations. For example, the worst-case residual error can be minimized via solving a min-max estimation problem as in Result <ref>.
Algorithm <ref> provides a robust estimation of $\alpha_j^*$ via solving the following min-max estimation problem:
α̂_j  =   min_α_j^*

max_ [[][c],c/c, Δ_x      Δy_j ] _F ≤δ_m    (K_xa + Δ_x K_aa^1/2)α_j^*  - ( y_j + Δy_j ) ,

where the upper-bound on the norm of perturbation matrix is given as
δ_m  =

√(1+  K_aa^1/2 ( K_xa^T K_xa   + τK_aa   +    β̂^-1 K_aa  )^-1  K_xa^T y_j  ^2)/ ( (τ+ β̂^-1)I + K_xa K_aa^-1K_xa^T   )^-1 y_j   ,

where $ \hat{\beta}^{-1}$ is the unique fixed point of $\mathcal{R}(\beta^{-1})$ to which the iterations (<ref>) converge.

The proof is provided in Appendix <ref>.

§ SECURE DISTRIBUTED DEEP LEARNING

§.§ Fuzzy Attributes and Classification Applications

[fuzzy attribute $\mathbf{A}_{\mathcal{P}_1}$] [fuzzy attribute $\mathbf{A}_{\mathcal{P}_2}$] [fuzzy attribute $\mathbf{A}_{\mathcal{P}_3}$]  [decision boundaries determined by the classifier (<ref>)]
A few examples of fuzzy attributes induced by deep autoencoders. The given set of 100 two-dimensional data samples, marked as `+' in the figures, has been used to train a wide conditionally deep membership-mapping autoencoder using Algorithm <ref> taking $L = 5$, $n = 2$, and possible $r_{max}$ values as $\{0.01,0.02,\cdots,0.1\}$. The membership function associated to a fuzzy attribute has been defined by (<ref>) taking $\nu = 2.001$.

A fuzzy attribute $\mathbf{A}_{\mathcal{P}}$, associated to a wide conditionally deep membership-mapping autoencoder $\mathcal{P}$ (that has been learned using dataset $\mathbf{Y} \subset \mathbb{R}^p$), can be defined on a universe of discourse $\mathbb{R}^p$ as
𝐀_𝒫     : =     { (y,μ_𝐀_𝒫(y))  |   y ∈ℝ^p }

where $\mu_{\mathbf{A}_{\mathcal{P}}}(y):\mathbb{R}^p \rightarrow [0,1]$ is a $p-$variate membership function such that $\mu_{\mathbf{A}_{\mathcal{P}}}(y)$ is interpreted as the degree to which a point $y \in \mathbb{R}^p$ matches to the attribute $\mathbf{A}_{\mathcal{P}}$. Without the loss of generality, the following Student-t type membership function can be defined to characterize $ \mathbf{A}_{\mathcal{P}} $:
μ_𝐀_𝒫(y)      =     (1 + 1/ν- 2   y - 𝒲𝒟(y;𝒫)^2 )^- ν+ p/2,

where $\widehat{\mathcal{WD}}(y;\mathcal{P})$ is through autoencoder $\mathcal{P}$ filtered output (<ref>). Similarly, one can define Gaussian type of membership function:
μ_𝐀_𝒫(y)      =     exp(- 1/2p  y - 𝒲𝒟(y;𝒫)^2  ).

A wide conditionally deep membership-mapping autoencoder induces a fuzzy attribute as defined in Definition <ref>. Fig. <ref> provides three different examples of fuzzy attributes induced by the deep autoencoders. As demonstrated through color-plots in Fig. <ref>, the defined fuzzy attribute (Definition <ref>) learns a representation of the data samples. This motivates to define a fuzzy classifier based on the following if-then rules:
,     ;

,     .

The class-label associated to a data point $y$ is predicted based on fuzzy rules (<ref>) as
𝒞(y;{ 𝒫_c}_c=1^C)     =       max_1 ≤c ≤C   μ_𝐀_𝒫_c(y).

The classifier (<ref>), $\mathcal{C}:\mathbb{R}^p \rightarrow \{1,2,\cdots,C \}$, assigns to an input vector the label of the class to which the data point has highest degree of matching. An example of the decision boundary determined by a three-class classifier is provided in Fig. <ref>.

§.§ Distributed Learning

We consider a scenario that data are distributed amongst different parties. Assume that there are $K$ different datasets, $\{\mathbf{Y}^1,\cdots,\mathbf{Y}^K \}$, owned locally by $K$ different parties. We consider the multi-class classification problem assuming that each local dataset, say $\mathbf{Y}^k$, can be partitioned into $C$ different classes, i.e.,
𝐘^k     =     { 𝐘^k_1,⋯, 𝐘^k_C }

where $\mathbf{Y}^k_c$ refers to the $c-$th class labelled data samples owned locally by the $k-$th party. Let $\mathcal{P}_c^k$ be the wide conditionally deep membership-mapping autoencoder learned from $\mathbf{Y}^k_c$ and $\mathbf{A}_{\mathcal{P}_c^k}$ be the corresponding fuzzy attribute. The fuzzy classifier (<ref>) can be extended for distributed setting as follows:
,     ;

,     .

The class-label associated to a data point $y$ is predicted based on fuzzy rules (<ref>) as
ĉ     =       max_1 ≤c ≤C   ( max_1 ≤k ≤K μ_ 𝐀_𝒫_c^k(y) )

=       min_1 ≤c ≤C   ( min_1 ≤k ≤K ( 1 - μ_ 𝐀_𝒫_c^k(y) ))

=       min_1 ≤c ≤C   ( min_1 ≤k ≤K  μ̅_ 𝐀_𝒫_c^k(y) ),

where $\bar{\mu}_{ \mathbf{A}_{\mathcal{P}_c^k}}(y) = 1 - \mu_{ \mathbf{A}_{\mathcal{P}_c^k}}(y)$.

§.§ Secure Homomorphic Evaluation of Global Classifier

[green,line width = 0.25mm](-3,-5)–(3,-5)–(3,3)–(-3,3)–cycle;
(0,-4.5) node[]$\begin{array}{c} \mbox{\scriptsize Party 1} \\ \mbox{\scriptsize \faMale\:\faFemale} \end{array}$;
(0,-3.25) node[rounded corners,draw=red,fill=red!5](n1)$\begin{array}{c}\mbox{local training data} \\ \mbox{to be protected} \\ \mbox{\small $\mathbf{Y}^1  =  \{ \mathbf{Y}^1_1,\cdots, \mathbf{Y}^1_C \}$} \end{array}$;
(0,-1.5) node[rounded corners,draw](n2) $\begin{array}{c}\mbox{local models} \\ \mbox{\small  $ \{ \mathcal{P}_1^1,\cdots, \mathcal{P}_C^1 \}$} \end{array}$;
[-latex,line width=0.2mm] (n1) to [out=90,in=-90] (n2);
(0,0.25) node[rounded corners,draw](n4) $\begin{array}{c}\mbox{inference} \\ \mbox{\small $\{\mathrm{cl}\_1,\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{\mathrm{cl}\_1}^1}}(y)\}$} \end{array}$;
[-latex,line width=0.2mm] (n2) to [out=90,in=-90] (n4);
[yellow,line width = 0.25mm](9,-5)–(15,-5)–(15,3)–(9,3)–cycle;
(11.75,-4.5) node[]$\begin{array}{c} \mbox{\scriptsize Party K} \\ \mbox{\scriptsize \faMale\:\faFemale} \end{array}$;
(12,0.25) node[rounded corners,draw](n11) $\begin{array}{c}\mbox{inference} \\ \mbox{\small $\{\mathrm{cl}\_K,\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{ \mathrm{cl}\_K}^K}}(y)\}$} \end{array}$;
(12,-1.5) node[rounded corners,draw](n14) $\begin{array}{c}\mbox{local models} \\ \mbox{\small $ \{ \mathcal{P}_1^K,\cdots, \mathcal{P}_C^K \}$} \end{array}$;
[-latex,line width=0.2mm] (n14) to [out=90,in=-90] (n11);
(12,-3.25) node[rounded corners,draw=red,fill= red!5](n15)$\begin{array}{c}\mbox{local training data} \\ \mbox{to be protected} \\ \mbox{\small $\mathbf{Y}^K  =  \{ \mathbf{Y}^K_1,\cdots, \mathbf{Y}^K_C \}$} \end{array}$;
[-latex,line width=0.2mm] (n15) to [out=90,in=-90] (n14);
[blue!40,line width = 0.25mm](4.25,-5)–(7.75,-5)–(7.75,2)–(4.25,2)–cycle;
(6,-4.5) node[]$\begin{array}{c} \mbox{\scriptsize User} \\ \mbox{\scriptsize \faUser} \end{array}$;
(6,-3.5) node[](n20) $\begin{array}{c}\mbox{input data} \\ \mbox{\small $y$} \end{array}$;
[thick,line width=0.2mm](n20) – (6,-2.5);
[thick,line width=0.2mm](3.5,-2.5) – (8.5,-2.5);
[-latex,thick,line width=0.2mm] (3.5,-2.5) to [out=90,in=0] (n4);
[-latex,thick,line width=0.2mm] (8.5,-2.5) to [out=90,in=180] (n11);
(0,2.1) node[rounded corners,draw, fill = gray!5](n7) $\begin{array}{c}\mbox{encryption {\scriptsize \faKey}} \\ \mbox{\small $\{\mathfrak{cp}_{\mathfrak{sk},N_b}(\mathrm{cl}\_1), \mathfrak{cp}_{\mathfrak{sk},N_b}(\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{\mathrm{cl}\_1}^1}}(y))\}$} \end{array}$;
[-latex,line width=0.2mm] (n4) to [out=90,in=-90] (n7);
(12,2.1) node[rounded corners,draw, fill = gray!5](n9) $\begin{array}{c}\mbox{encryption {\scriptsize \faKey}} \\ \mbox{\small $\{\mathfrak{cp}_{\mathfrak{sk},N_b}(\mathrm{cl}\_K), \mathfrak{cp}_{\mathfrak{sk},N_b}(\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{\mathrm{cl}\_K}^K}}(y))\}$} \end{array}$;
[-latex,line width=0.2mm] (n11) to [out=90,in=-90] (n9);
draw = gray,
fill = cyan!5,
minimum width = 11cm,
minimum height = 5cm,
cloud puffs = 18] (c) at (6,5) ;
(6,4.75) node[rounded corners,draw, fill=gray!5](n8) $\begin{array}{c}\mbox{global model inference} \\ \mbox{\small $\left\{\mathfrak{cp}_{\mathfrak{sk},N_b}(\mathrm{cl}\_k), \mathfrak{ct}_{\mathfrak{sk}}(\delta(\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{\mathrm{cl}\_k}^k}}(y), \bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{\mathrm{cl}\_{k^*}}^{k^*}}}(y) ))\: | \: k = 1,\cdots,K\right \}$} \end{array}$;
[-latex,line width=0.2mm,cyan] (n7) to [out=90,in=180] (n8);
(6,6.25) node[]$\begin{array}{c} \mbox{homomorphic evaluation} \\ \mbox{of global model} \end{array}$  ;
[-latex,line width=0.2mm,cyan] (n9) to [out=90,in=0] (n8);
[orange,line width = 0.25mm] (4.5,-0.75)–(7.5,-0.75)–(7.5,1.5)–(4.5,1.5)–cycle;
(6,0.5) node[rounded corners,draw, fill = gray!5](n17) $\begin{array}{c}\mbox{output data} \\ \mbox{\small $\hat{c}$} \end{array}$;
(6,-0.4) node[]decryption ;
[-latex,line width=0.2mm,cyan] (n8) to [out=-90,in=90] (n17);

A practical method for secure distributed deep learning based on membership-mappings and fully homomorphic encryption.

For a secure distributed learning based on fully homomorphic encryption, a practical approach is suggested based on the observation that instead of encrypting higher-dimensional vectors, it is sufficient to encrypt the scalar values $\{\bar{\mu}_{ \mathbf{A}_{\mathcal{P}_c^k}}(y) \; | \; c = 1,\cdots,C; \: k = 1,\cdots,K\}$ for the homomorphic evaluation of classifier (<ref>) (i.e. evaluation of (<ref>)) in an efficient manner. That is, the homomorphic evaluation of the “$\arg \; \min$” function over $K \times C$ encrypted values is required. More efficiently, the homomorphic evaluation of the “$\arg \; \min$” function over $K$ encrypted values is sufficient for the homomorphic evaluation of the classifier. For this, define $\mathrm{cl}\_k$ as the class-label predicted by the $k-$th local classifier, i.e.,
cl_k     =        𝒞(y;{ 𝒫_c^k}_c=1^C)

where $ \mathcal{C}(\cdot)$ is defined as in (<ref>). Now, (<ref>) can be alternatively expressed as
ĉ     =     cl_k^*

k^*     =       min_1 ≤k ≤K    μ̅_𝐀_ 𝒫_ cl_k^k(y).

Let $\delta(\mathfrak{m}_1,\mathfrak{m}_2)$ be the Kronecker delta function of two variables $\mathfrak{m}_1,\mathfrak{m}_2 \in [0,1]$ defined as
δ(𝔪_1,𝔪_2)     =     {
1     if $\mathfrak{m}_1=\mathfrak{m}_2$,

0     if $\mathfrak{m}_1 \neq \mathfrak{m}_2$.

It follows from (<ref>-<ref>) that
ĉ     =     ∑_k=1^K cl_k   δ(μ̅_𝐀_ 𝒫_ cl_k^k(y),μ̅_𝐀_ 𝒫_ cl_k^*^k^*(y)).

A practical method using TFHE scheme <cit.> is provided for secure distributed deep learning. For this, define the followings:

* For a given positive integer $N_b \in \mathbb{Z}_{>0}$, let $\mathfrak{pt}_{N_b}:[0,1]\rightarrow \{0,1,\cdots,2^{N_b}-1 \}$ be a function defined as
𝔭𝔱_N_b(𝔪)     :=     ⌈(2^N_b-1) 𝔪 ⌉,   𝔪 ∈[0,1].

In our setting, $\mathfrak{pt}_{N_b}( \mathfrak{m})$ is the plaintext that encodes a message $ \mathfrak{m} \in [0,1]$ as unsigned $N_b-$bit integer.

* Let $\mathrm{BitDec}_{N_b}: \{0,1,\cdots,2^{N_b}-1 \} \rightarrow \{0,1\}^{N_b}$ be the binary representation of a $N_b-$bit unsigned integer. That is,
(𝔟𝔱_1( 𝔪),⋯,𝔟𝔱_N_b( 𝔪))      =     BitDec_N_b(𝔭𝔱_N_b( 𝔪)),
where $\mathfrak{bt}_k( \mathfrak{m}) \in \{0,1\}$ for all $k \in \{1,2,\cdots,N_b\}$.

* Let $N_c$ be the ciphertext dimension set for a given value of security bits, say 128 bits security.

* Let $\mathfrak{sk} \in \{0,1\}^{N_c}$ be a secret key generated for TFHE encryption.

* Let $\mathfrak{ct}_{\mathfrak{sk}}(\mathfrak{bt}) \in  \mathbb{T}^{N_c+1}$, where $\mathbb{T} = \mathbb{R}/\mathbb{Z}$, be the TFHE encryption of a bit $\mathfrak{bt} \in \{0,1 \}$, i.e., $\mathfrak{ct}_{\mathfrak{sk}}(\mathfrak{bt}) = \mathrm{TFHE.Enc}(\mathfrak{bt};\mathfrak{sk})$.

* Let $\mathfrak{cp}_{\mathfrak{sk},N_b}:[0,1] \rightarrow \mathbb{T}^{N_b (N_c+1)}$ be a function defined as
𝔠𝔭_𝔰𝔨,N_b(𝔪)     :=     (𝔠𝔱_𝔰𝔨(𝔟𝔱_1( 𝔪)),⋯,𝔠𝔱_𝔰𝔨(𝔟𝔱_N_b( 𝔪)))
where $\mathfrak{ct}_{\mathfrak{sk}}\left(\mathfrak{bt}_k( \mathfrak{m})\right)$ is the TFHE encryption of bit $\mathfrak{bt}_k( \mathfrak{m})$. Thus, $\mathfrak{cp}_{\mathfrak{sk},N_b}(\mathfrak{m})$ homomorphically encrypts the message $\mathfrak{m} \in [0,1]$ with $N_b$-bit precision.

Our approach to homomorphically evaluate the global classifier (<ref>) is shown in Fig. <ref>. The approach consists of following steps:

* A pair of secret and cloud keys is generated. The secret key is meant for encryption and decryption. The cloud key is exported to the cloud, and allows to operate over encrypted data.

* For a given input $y$, the outputs of local classifiers $\{\mathfrak{cp}_{\mathfrak{sk},N_b}(\mathrm{cl}\_k), \mathfrak{cp}_{\mathfrak{sk},N_b}(\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{\mathrm{cl}\_k}^k}}(y))\; | \; k = 1,\cdots,K\}$ are sent to the cloud for performing secure homomorphic computations.

* The values $\{  \delta (\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{ \mathrm{cl}\_k}^k}}(y),\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{ \mathrm{cl}\_{k^*}}^{k^*}}}(y))  \; |  \: k = 1,\cdots,K\}$ are homomorphically evaluated in the cloud from the encrypted data (sent by different parties) using a boolean circuit composed of bootstrapped binary gates.

* The encrypted output of the global model $\{\mathfrak{cp}_{\mathfrak{sk},N_b}(\mathrm{cl}\_k), \mathfrak{ct}_{\mathfrak{sk}}(\delta(\bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{\mathrm{cl}\_k}^k}}(y), \bar{\mu}_{\mathbf{A}_{ \mathcal{P}_{\mathrm{cl}\_{k^*}}^{k^*}}}(y) ))\: | \: k = 1,\cdots,K\}$ is sent to the user.

* The class-label associated to a data point $y$ is predicted using (<ref>) after decrypting the data provided by the cloud.

§ EXPERIMENTS

The proposed method was implemented using MATLAB R2017b and TFHE C/C++ library <cit.> on a MacBook Pro machine with a 2.2 GHz Intel Core i7 processor and 16 GB of memory. The previous works <cit.> have already verified the competitive performance of membership-mappings in classification applications. In this study, our focus is to verify the application potential of the proposed approach and to compare the method with the alternative differentially private distributed deep learning approach <cit.>.

§.§ Details

Under a distributed deep learning scenario, the local models are developed using Algorithm <ref> taking $L = 5$, $n = 20$, and possible $r_{max}$ values array as $\{0.5\}$. Targeting 128-bits of security, TFHE library is used to homomorphically evaluate the global classifier with the precision of 16-bits and also 8-bits. For a comparison, differentially private local models, developed using Algorithm <ref> on the noisy training data obtained via optimal noise adding mechanism <cit.> (taking adjacency parameter $d = 1$, failure probability $\delta = 1\mathrm{e}{-5}$, and privacy-loss bound $\epsilon = 0.1$ and also $\epsilon = 1$), are combined through fuzzy rules (<ref>) without any encryption. The test data accuracy and average computational time required for secure homomorphic computations in the cloud (for computing the encrypted global output for a given input) are considered as performance indices.

§.§ MNIST Dataset

The first experiment is on the widely used MNIST digits dataset containing $28 \times 28$ sized images divided into training set of 60000 images and testing set of 10000 images. The images' pixel values were divided by 255 to normalize the values in the range from $0$ to $1$. The $28 \times 28$ normalized values of each image were flattened to an equivalent $784-$dimensional data vector. A two-party scenario is considered such that Party-A owns all the training images of odd digits while Party-B owns the rest training images of even digits.

Experimental Results on MNIST Dataset

Method     Accuracy     Time

$\begin{array}{c} \mbox{Proposed (8-bits precision of} \\ \mbox{homomorphic computations)} \end{array}$     0.9762     3.2251 s

$\begin{array}{c} \mbox{Proposed (16-bits precision of} \\ \mbox{homomorphic computations)} \end{array}$     0.9872     4.9785 s

$\begin{array}{c} \mbox{differentially private} \\ \mbox{distributed deep learning with} \\ \mbox{privacy-loss bound $\epsilon = 0.1$ \cite{KUMAR202187}} \end{array}$     0.0892     n/a

$\begin{array}{c} \mbox{differentially private} \\ \mbox{distributed deep learning with} \\ \mbox{privacy-loss bound $\epsilon = 1$ \cite{KUMAR202187}} \end{array} $     0.8994     n/a

Table <ref> reports the experimental results.

§.§ Freiburg Groceries Dataset

The second experiment is on “Freiburg Groceries Dataset” considered previously <cit.> for privacy-preserving distributed learning experiments. The dataset contains 4947 labeled images of grocery products categorized into 25 different classes. A feature vector is created from each image by extracting features from “AlexNet” and “VGG-16” networks which are pre-trained Convolutional Neural Networks. The activations of the fully connected layer “fc6” in AlexNet constitute a $4096-$dimensional feature vector. Similarly, the activations of the fully connected layer “fc6” in VGG-16 constitute another $4096-$dimensional feature vector. The features extracted by both networks are joined together to form a $8192-$dimensional vector. The feature vectors are normalized to have zero-mean and unity-variance along each dimension. The set of normalized feature vectors is split into a training set containing around 80% of data points and a testing set containing remaining data points. A three-party scenario is created such that Party-A owns all the training data of first ten grocery categories, Party-B owns all the training data of second ten grocery categories, and Party-C owns all the training data of rest five grocery categories.

Experimental Results on Freiburg Groceries Dataset

Method     Accuracy     Time

$\begin{array}{c} \mbox{Proposed (8-bits precision of} \\ \mbox{homomorphic computations)} \end{array}$     0.8861     4.9534 s

$\begin{array}{c} \mbox{Proposed (16-bits precision of} \\ \mbox{homomorphic computations)} \end{array}$     0.8880      7.9240 s

$\begin{array}{c} \mbox{differentially private} \\ \mbox{distributed deep learning with} \\ \mbox{privacy-loss bound $\epsilon = 0.1$ \cite{KUMAR202187}} \end{array}$     0.1356     n/a

$\begin{array}{c} \mbox{differentially private} \\ \mbox{distributed deep learning with} \\ \mbox{privacy-loss bound $\epsilon = 1$ \cite{KUMAR202187}} \end{array}$     0.8261     n/a

Table <ref> reports the experimental results.

§.§ A Biomedical Application

As an application example, the mental stress detection problem is considered. The dataset from <cit.>, consisting of heart rate interval measurements of different subjects, is considered for the study of individual stress detection problem. The problem is concerned with the detection of stress on an individual based on the analysis of recorded sequence of R-R intervals, $\{RR^i\}_i$. The R-R data vector at $i-$th time-index, $y^i$, is defined as $y^i  =  \left[\begin{IEEEeqnarraybox*}[][c]{,c/c/c/c,} RR^i & RR^{i-1} & \cdots & RR^{i-d}\end{IEEEeqnarraybox*} \right]^T$. That is, the current interval and history of previous $d$ intervals constitute the data vector. Assuming an average heartbeat of 72 beats per minute, $d$ is chosen as equal to $72 \times 3 = 216$ so that R-R data vector consists of on an average 3-minutes long R-R intervals sequence. A dataset, say $\{y^i\}_{i}$, is built via 1) preprocessing the R-R interval sequence $\{RR^i\}_i$ with an impulse rejection filter for artifacts detection, and 2) excluding the R-R data vectors containing artifacts from the dataset. The dataset contains the stress-score on a scale from 0 to 100. A label of either “no-stress” or “under-stress” is assigned to each $y^i$ based on the stress-score. Thus, we have a binary classification problem. A two-party collaborative learning scenario is considered where a randomly chosen subject is considered as Party-A. While keeping Party-A fixed, the distributed learning experiments are performed independently on every other subject being considered as Party-B. For each subject, 50% of the data samples serve as training data while remaining as test data. The subjects, with data containing both the classes and at least 60 samples, are considered for experimentation. There are in total 48 such subjects.

Experimental Results on Heart Rate Variability Dataset

Method     Accuracy     Time

$\begin{array}{c} \mbox{Proposed (8-bits precision of} \\ \mbox{homomorphic computations)} \end{array}$     0.8358      3.1571 s

$\begin{array}{c} \mbox{Proposed (16-bits precision of} \\ \mbox{homomorphic computations)} \end{array}$      0.9580     4.8610 s

$\begin{array}{c} \mbox{differentially private} \\ \mbox{distributed deep learning with} \\ \mbox{privacy-loss bound $\epsilon = 0.1$ \cite{KUMAR202187}} \end{array}$      0.5123     n/a

$\begin{array}{c} \mbox{differentially private} \\ \mbox{distributed deep learning with} \\ \mbox{privacy-loss bound $\epsilon = 1$ \cite{KUMAR202187}} \end{array}$     0.6873     n/a

The experimental results, averaged over 48 independent experiments, are reported in Table <ref>.

§.§ Scalability

The computational time required for the homomorphic evaluation of global model depends on the number of parties (i.e. $K$) participating in collaborative learning. Therefore, experiments are performed to study the computational time as the number of parties is varied from $K = 2$ to $K = 100$.

Computational time vs. number of parties.

Fig. <ref> plots the computational time required for secure homomorphic computations on a MacBook Pro machine with a 2.2 GHz Intel Core i7 processor and 16 GB of memory without using parallel computing. The ratio of computational time to the number of parties is observed to be approximately equal to 3 s for 16-bit precision and 1.75 s for 8-bit precision.

§.§ Evaluation of the Proposed Approach

To evaluate the gains achieved in accuracy and computational time as a result of the proposed approach to secure distributed deep learning (summarized in Fig. <ref>), a comparison is made with a variant of the TFHE fully homomorphic encryption scheme <cit.>. The study in <cit.> reports the experiments on MNIST dataset for evaluating the neural networks with different depths (referred to as NN-20, NN-50, and NN-100) over TFHE fully homomorphically encrypted data. The results of <cit.> are compared with the proposed method (considering two-party scenario with Party-A owning odd-digit images and Party-B owning even-digit images) in Table <ref>.

Gains achieved by the proposed approach on MNIST dataset

Method     Security     Machine     Accuracy     Time

Proposed     128-bits     $\begin{array}{c} \mbox{2.2 GHz} \\ \mbox{Intel Core i7} \end{array}$     0.987     4.98 s

NN-20 <cit.>     128-bits     $\begin{array}{c} \mbox{2.6 GHz} \\ \mbox{Intel Core i7} \end{array}$     0.971       115.52 s

NN-50 <cit.>     128-bits     $\begin{array}{c} \mbox{2.6 GHz} \\ \mbox{Intel Core i7} \end{array}$     0.947       233.55 s

NN-100 <cit.>     128-bits     $\begin{array}{c} \mbox{2.6 GHz} \\ \mbox{Intel Core i7} \end{array}$     0.830       481.61 s

§.§ Main Results

Following inferences are made from the results of aforementioned experiments.

* As expected, the proposed approach leads to better accuracy (observed in Tables <ref>, <ref>, <ref>) than the differentially private approach <cit.>, since differential privacy requires contaminating data with noise to preserve privacy.

* The proposed membership-mappings based approach is capable of handling the large computational overhead issue of fully homomorphic encryption, since the computational time on a MacBook Pro machine with a 2.2 GHz Intel Core i7 processor and 16 GB of memory (as reported in Tables <ref>, <ref>, <ref>, <ref>) is practical.

* A linear increase of the computational time with increasing number of parties, as observed in Fig. <ref>, indicates that the proposed approach is scalable using parallel computing.

* Remarkably, the computational time required for secure homomorphic evaluation of the global model in the cloud is independent of the dimension of the input data, and thus the approach is practical. The computational time depends only on the number of parties and the chosen precision. The ratio of computational time to the number of parties as approximately equal to 3 s for 16-bit precision verifies the application potential.

§ CONCLUDING REMARKS

This study has outlined a membership-mappings based approach to secure distributed deep learning. The crux of our methodology lies in defining fuzzy attributes (which are induced by globally convergent and robust variational membership-mappings based local deep models) allowing to combine local models by means of a rule-based fuzzy system, thus facilitating the homomorphic evaluation of the global model efficiently. The feature, that the computational time for secure homomorphic evaluation of the global model in the cloud is independent of the dimension of input data, adds to the practicality of the approach. The experimental results verify that the proposed method, while preserving the privacy in a distributed learning scenario using fully homomorphic encryption, remains accurate, practical, and scalable.

§ MEMBERSHIP-MAPPING OUTPUT ESTIMATION

Using (<ref>) and (<ref>), we have
< (f_j)_i>_μ_f_j;u_j       =     (K_xa K_aa^-1  u_j )_i

=     G(x^i) K_aa^-1  u_j.

ℱ_j(x^i)  =      G(x^i)K_aa^-1 < u_j >_μ^*_u_j.

Using (<ref>) and (<ref>) in (<ref>), we have
ℱ_j(x^i)      =   β(G(x^i) ) K_aa^-1 K̂_u_j K_aa^-1K_xa^T y_j .

Substituting $\hat{K}_{\mathrm{u}_j}$ from (<ref>) in (<ref>), we get (<ref>).

§ PROOF OF RESULT <REF>

The proof is split into three parts.

Part 1
Since $K_{\mathrm{a}\mathrm{a}} > 0$, there exists the unique square root, $K_{\mathrm{a}\mathrm{a}}^{1/2} > 0$. Thus,
S V^T  K_aa^-1 VS     =     (K_aa^-1/2VS )^T( K_aa^-1/2VS )

>     0.

Since $\tau > 0$ and $\beta > 0$,
min_eigen(I + 1/(τ+ β^-1) S V^T  K_aa^-1 VS )     >     1

where “$\text{min\_eigen}(\cdot)$” denotes the minimum eigenvalue. Thus,
max_eigen((I + 1/(τ+ β^-1) S V^T  K_aa^-1 VS )^-2 )     <     1
where “$\text{max\_eigen}(\cdot)$” denotes the maximum eigenvalue. As a result of (<ref>),
ℛ(β^-1)      <      1/pN∑_j=1^p(  b_j^1 ^2 +  b_j^2 ^2).

As $U$ is orthogonal, it follows from (<ref>) that $ \|  b_j^1 \|^2 +\|  b_j^2 \|^2 = \| \mathrm{y}_j  \|^2$, and thus
ℛ(β^-1)      <      β^-1|_up.

It follows immediately from (<ref>) that $\mathcal{R}(\beta^{-1})   >   \beta^{-1}|_{low}$. Hence, $\mathcal{R}(\beta^{-1}) \in \left(\beta^{-1}|_{low}, \beta^{-1}|_{up}  \right)$.

Part 2
The derivative of $\mathcal{R}$ w.r.t. $\beta^{-1}$ is given as
ℛ(β^-1)/β^-1  =

2/pN∑_j=1^p { (τ+ β^-1)(b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2  b_j^1 .

.  -     (τ+ β^-1)^2 (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-3  b_j^1 }.

(τ+ β^-1)^2 (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-3  b_j^1

≤    (τ+ β^-1)^2  ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-1 _2

×   (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2 b_j^1

=     (τ+ β^-1)  ( I + 1/(τ+ β^-1) S V^T  K_aa^-1 VS   )^-1 _2

×   (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2 b_j^1

=      (τ+ β^-1) (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2 b_j^1/min_sing( I + 1/(τ+ β^-1) S V^T  K_aa^-1 VS )
where “$\text{min\_sing}(\cdot)$” denotes the minimum singular value. Observing that $\tau > 0$, $\beta^{-1} > 0$, and $S V^T  K_{\mathrm{a}\mathrm{a}}^{-1} VS > 0$ (i.e. (<ref>)), we have
min_sing( I + 1/(τ+ β^-1) S V^T  K_aa^-1 VS )

=     1 + min_sing(1/(τ+ β^-1) S V^T  K_aa^-1 VS )

>     1.

Combining (<ref>) and (<ref>), we have
(τ+ β^-1)^2 (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-3  b_j^1

<      (τ+ β^-1)(b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2 b_j^1.

Using (<ref>) in (<ref>), we get
ℛ(β^-1)/β^-1      >     0.

Observing that $\tau > 0$, $\beta^{-1} > 0$, and $S V^T  K_{\mathrm{a}\mathrm{a}}^{-1} VS > 0$ (i.e. (<ref>)), we also have
(τ+ β^-1)^2 (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-3  b_j^1     >     0.
Using (<ref>) in (<ref>), we get
ℛ(β^-1)/β^-1 <

2/pN (τ+ β^-1) ∑_j=1^p (b_j^1)^T ( (τ+ β^-1)I + S V^T  K_aa^-1 VS   )^-2  b_j^1.

Inequality (<ref>), using (<ref>), can be expressed as
ℛ(β^-1)/β^-1      <     2/pN pNℛ(β^-1)-∑_j=1^p    b_j^2 ^2/τ+β^-1

<     2/pN ∑_j=1^p    b_j^1 ^2/τ+β^-1,

where (<ref>) follows using (<ref>). Inequalities (<ref>) and (<ref>) lead to (<ref>).

Part 3
Introduce $h(\beta^{-1})  =  \mathcal{R}(\beta^{-1}) - \beta^{-1}$, and observe that $h(\beta^{-1}|_{low}) > 0$ and $h(\beta^{-1}|_{up}) < 0$. By the intermediate value theorem, there is a $\hat{\beta}^{-1} \in \left(\beta^{-1}|_{low}, \beta^{-1}|_{up}  \right)$ such that $h(\hat{\beta}^{-1}) = 0$, i.e., $\hat{\beta}^{-1}  =\mathcal{R}(\hat{\beta}^{-1})$. Thus, $\hat{\beta}^{-1}$ is a fixed point of $\mathcal{R}(\beta^{-1})$.

§ PROOF OF RESULT <REF>

As $\beta^{-1}|_{it} > \beta^{-1}|_{low}$, it follows from (<ref>) that
τ+ β^-1|_it     >      2/pN∑_j=1^p   y_j  ^2,

and thus
ℛ(β^-1|_it)/β^-1      <     ∑_j=1^p    b_j^1 ^2/∑_j=1^p   y_j  ^2
=     ∑_j=1^p    b_j^1 ^2/∑_j=1^p (  b_j^1 ^2 +    b_j^2 ^2).
That is, there exists a constant $k$ such that
0     <     ℛ(β^-1|_it)/β^-1      ≤    k     < 1,  ∀it ∈{0,1,2,⋯}.

Let $\hat{\beta}^{-1}$ be a fixed point of $\mathcal{R}(\beta^{-1})$. Now, consider
| β^-1|_it -  β̂^-1 |     =      | ℛ(β^-1|_it-1) - ℛ(β̂^-1) |

≤    k | β^-1|_it-1 -  β̂^-1 |

≤    k^it | β^-1|_0 - β̂^-1 |,

that leads to
lim_it →∞ | β^-1|_it -  β̂^-1 |     ≤    lim_it →∞ k^it | β^-1|_0 -  β̂^-1 |     =     0.
The uniqueness of the fixed point can be seen via assuming by contradiction that there exists another fixed point, say $\tilde{\beta}^{-1}$. Now consider
| β̃^-1 -  β̂^-1 |     =     |ℛ(β̃^-1) -  ℛ(β̂^-1) |     ≤     k | β̃^-1 -  β̂^-1 |

<     | β̃^-1 -  β̂^-1 |.
This implies that $\tilde{\beta}^{-1}  =  \hat{\beta}^{-1} $. Hence, the result follows.

§ PROOF OF RESULT <REF>

According to the triangle inequality,
(K_xa + Δ_x K_aa^1/2)α_j^* - ( y_j + Δy_j )

≤       K_xa α_j^* - y_j   +  Δ_x K_aa^1/2 α_j^* -  Δy_j

≤     K_xa α_j^* - y_j   +  [[][c],c/c, Δ_x      Δy_j ] _F  [[][c],c, K_aa^1/2 α_j^*
-1  ]

≤     K_xa α_j^* - y_j   + δ_m √(1+ (α_j^*)^T  K_aa α_j^*),

and hence
α̂_j     =     ( K_xa^T K_xa + δ_m   K_xa α̂_j - y_j  /√(1+ α̂_j^TK_aa α̂_j )K_aa )^-1 K_xa^T y_j.
Algorithm <ref> estimates $\alpha_j(\hat{\beta}^{-1})$ using (<ref>). It  can be seen using (<ref>) that
(τ+ β̂^-1) √(1+(α_j(β̂^-1))^TK_aaα_j(β̂^-1) )/   K_xa α_j(β̂^-1) - y_j         =     δ_m.

As a result of (<ref>), it follows from (<ref>) that
α_j(β̂^-1)  =  ( K_xa^T K_xa .

. +  δ_m    K_xa α_j(β̂^-1) - y_j   /√(1+(α_j(β̂^-1))^TK_aaα_j(β̂^-1) )   K_aa )^-1  K_xa^T y_j.
As equalities (<ref>) and (<ref>) are identical, $\alpha_j(\hat{\beta}^{-1})$ (which is the solution of (<ref>)) must be equal to $\hat{\alpha}_j$ (which is the solution of (<ref>)), i.e.,
α_j(β̂^-1)     =     α̂_j.

Hence, Algorithm <ref> solves the min-max problem (<ref>).