## 4 Proposed Method
### Lower Bound of Log-Posterior Distribution  
The inducing method [47] uses the property of the conditional distribution of a Gaussian and replaces the matrix inversion of a large matrix with that of a small matrix. Specifically, we define a few inducing points \(\mathbf{u}\in\mathbb{R}^{M}\) and their positions \(\mathbf{Z}=[\mathbf{z}_{1},\mathbf{z}_{2},\ldots,\mathbf{z}_{M}]^{\top}\in \mathbb{R}^{M\times Q}\) (\(M\) being the number of inducing points) in the latent space and evaluate the log-likelihood function after conditioning by these inducing positions \(\mathbf{Z}\) as \(\log p(\mathbf{y}_{,d}|\mathbf{X},\mathbf{Z})\). First, the joint distribution of \(\mathbf{y}_{,d}\), \(\mathbf{f}\), and \(\mathbf{u}\) is given as follows:  
\[p(\mathbf{y}_{,d},\mathbf{f},\mathbf{u}|\mathbf{X},\mathbf{Z})=p(\mathbf{y}_{,d}|\mathbf{f})p(\mathbf{f}|\mathbf{u},\mathbf{X},\mathbf{Z})p(\mathbf{u}| \mathbf{Z}). \tag{20}\]  
Each distribution is given as the following equation:  
\[p(\mathbf{y}_{,d}|\mathbf{f}) =\mathcal{N}(\mathbf{y}_{,d}|\mathbf{f},\beta^{-1}\mathbf{I}), \tag{21}\] \[p(\mathbf{f}|\mathbf{u},\mathbf{X},\mathbf{Z}) =\mathcal{N}(\mathbf{f}|\mathbf{K}_{NM}\mathbf{K}_{MM}\mathbf{u}, \mathbf{K}_{NN}-\mathbf{K}_{NM}\mathbf{K}_{MM}^{-1}\mathbf{K}_{MN}),\] (22) \[p(\mathbf{u}|\mathbf{Z}) =\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{MM}), \tag{23}\]  
where \(\mathbf{K}_{MN}=\mathbf{K}_{NM}^{\top}\in\mathbb{R}^{M\times N}\) denotes a covariance matrix whose \(ij\)-th entry is \(k_{SM}(\mathbf{x}_{i},\mathbf{z}_{j})\) and \(\mathbf{K}_{MM}\in\mathbb{R}^{M\times M}\) is a gram matrix whose \(ij\)-th entry denotes \(k_{SM}(\mathbf{z}_{i},\mathbf{z}_{j})\). For the evaluation of \(\log p(\mathbf{y}_{,d}|\mathbf{X},\mathbf{Z})\), we need to marginalize the Gaussian process prior \(\mathbf{f}\) and inducingpoints \(\mathbf{u}\). The variational method [47] introduces the Jensen inequality to evaluate the log-likelihood as follows:  
\[\log p(\mathbf{y}_{:d}|\mathbf{X},\mathbf{Z}) =\log\int\biggl{\{}\int p(\mathbf{y}_{:d}|\mathbf{f})p(\mathbf{f}| \mathbf{u},\mathbf{X},\mathbf{Z})d\mathbf{f}\biggr{\}}p(\mathbf{u}|\mathbf{Z}) d\mathbf{u}\] \[\geq\log\int\exp\bigg{\{}\int\log p(\mathbf{y}_{:d}|\mathbf{f})p( \mathbf{f}|\mathbf{u},\mathbf{X},\mathbf{Z})d\mathbf{f}\biggr{\}}p(\mathbf{u} |\mathbf{Z})d\mathbf{u}\] \[\triangleq\mathcal{L}_{d}. \tag{24}\]  
We evaluate the lower bound \(\mathcal{L}_{d}\) of the log-posterior in Equation (24) instead of the original log-likelihood. \(\mathcal{L}_{d}\) can be computed in the closed-form expression as follows:  
\[\mathcal{L}_{d}=-\frac{N}{2}\log 2\pi-\frac{1}{2}\log|\mathbf{Q}_{NN}+\beta^{- 1}\mathbf{I}|-\frac{1}{2}\mathbf{y}_{:d}^{\top}(\mathbf{Q}_{NN}+\beta^{-1} \mathbf{I})^{-1}\mathbf{y}_{:d}-\frac{\beta}{2}\text{tr}(\mathbf{K}_{NN}- \mathbf{Q}_{NN}), \tag{25}\]  
where \(\mathbf{Q}_{NN}=\mathbf{K}_{NM}\mathbf{K}_{MM}^{-1}\mathbf{K}_{MN}\). By comparing Equation (4) with Equation (25), the \(N\times N\) matrix \(\mathbf{K}_{NN}\) is replaced by the matrix \(\mathbf{Q}_{NN}\), which is considered a low-rank approximated form of the original matrix \(\mathbf{K}_{NN}\). Therefore, by applying the matrix formula to Equation (25), we can evaluate the lower bound of the likelihood in \(O(NM^{2})\) time complexity more efficiently than the likelihood function in Equation (4). By substituting (25) into (17) and neglecting the constants, we obtain the following objective function \(\mathcal{L}\):  
\[\mathcal{L}=-\frac{D}{2}\log|\mathbf{Q}_{NN}+\beta^{-1}\mathbf{ I}|-\frac{1}{2}\text{tr}\Bigl{[} (\mathbf{Q}_{NN}+\beta^{-1}\mathbf{I})^{-1}\mathbf{Y}\mathbf{Y}^{\top} \Bigr{]}\\ -\frac{\beta D}{2}\text{tr}(\mathbf{K}_{NN}-\mathbf{Q}_{NN})- \frac{\alpha Q}{2}\text{tr}\Bigl{(}\mathbf{L}\mathbf{X}\mathbf{X}^{\top} \Bigr{)}. \tag{26}\]  
We maximize the lower bound of the likelihood function in Equation (26) and explain the methodology of the optimization in the following section.