# Enhancing Diffusion-based Point Cloud Generation with Smoothness Constraint

Yukun Li
Department of Computer Science
Tufts University
Medford, MA 02155
<EMAIL_ADDRESS>
&Liping Liu
Department of Computer Science
Tufts University
Medford, MA 02155
<EMAIL_ADDRESS>

###### Abstract

Diffusion models have been popular for point cloud generation tasks. Existing
works utilize the forward diffusion process to convert the original point
distribution into a noise distribution and then learn the reverse diffusion
process to recover the point distribution from the noise distribution.
However, the reverse diffusion process can produce samples with non-smooth
points on the surface because of the ignorance of the point cloud geometric
properties. We propose alleviating the problem by incorporating the local
smoothness constraint into the diffusion framework for point cloud generation.
Experiments demonstrate the proposed model can generate realistic shapes and
smoother point clouds, outperforming multiple state-of-the-art methods.

##  1 Introduction

As a widely-used 3D representation, point clouds have recently attracted
widespread attention due to their compactness and proximity to raw sensory
data [1, 2, 3, 4]. Several sophisticated algorithms have been proposed for
analyzing point clouds in robotics and autonomous driving [5, 6, 7, 8, 9].
However, acquiring raw data in real applications demands significant labor
effort and incurs time costs. Moreover, the obtained point clouds are
frequently imperfect; they may be sparse and partial due to factors such as
distances, occlusions, reflections, and the limitations of device resolution
and angles [10, 11]. Recently, diffusion-based generative models have been
effective in generating 2D images and pose a promising direction for
generating 3D point clouds. However, point cloud generation poses more
significant challenges than image generation since points in a point cloud
have no order and are arranged irregularly. In addition, finding a way to
ensure smoothness on the surface of the point cloud is also a challenging
problem.

Several algorithms have been proposed for point cloud generation. For
instance, probabilistic methods [12, 13] take the generation of point clouds
as a probabilistic inference problem. Building upon the structure of
Variational Autoencoders (VAE), these methods aim to learn the distribution of
shapes. The process involves initially learning the distribution of a shape
and subsequently learning the conditional distribution of points generated
based on this shape. But the Pointflow [12] needs to normalize the learned
distribution, which is time-consuming. ShapeGF [14] must require two stages of
training. Some methods are based on the autoregressive model [15], but these
methods need to specify an order of points during training. [16] proposed a
point cloud GAN model with a hierarchical sampling to generate new point cloud
objects. Other GAN-based methods like [17, 18, 19] are challenging to train
and unstable, primarily due to the use of the adversarial loss function during
the training process.

There are also several diffusion-based generative models [13, 20], which treat
points as particles in a thermodynamic system that diffuse from the original
distribution to the noise distribution. The generation process can move the
noisy point locations toward the original point cloud surface. The training
aims to learn a score model to approximate the gradient of log-likelihood of
the point cloud distribution, which is equivalent to learning the moving
direction of the particles (points) during diffusion. However, both algorithms
are discrete generative processes and thus cannot consider the shape of
geometric information when learning point-moving directions. Compared with
diffusion-based models for image data that learn the data distribution in the
feature space, these models face a critical challenge in learning the 3D point
location distribution. This paper aims to incorporate the local geometric
constraints into a diffusion-based point cloud generation process to enforce
surface smoothness. Our approach views the diffusion probabilistic model from
a Bayesian perspective to consider the prior constraint of the point cloud.
Experiments demonstrate that our method is not only capable of generating
realistic shapes but also generating smoother point clouds, outperforming
multiple state-of-the-art methods.

##  2 Background

###  2.1 Tweedie’s Formula

Consider the noisy sample distribution
$p(\hat{\mathbf{X}}|\mathbf{X})\sim\mathcal{N}(\mathbold{\mu},\sigma^{2}\mathbf{I})$,
the Tweedie’s formula [21, 22] tells us the optimum denoised sample can be
achieved by computing the posterior expectation:

$\displaystyle\mathbb{E}[\mathbf{X}|{\hat{\mathbf{X}}}]=\hat{\mathbf{X}}+\sigma^{2}\nabla_{\hat{\mathbf{X}}}\log
p(\hat{\mathbf{X}})$ (1)

Thus, for a diffusion model in which the forward step is modeled as
$p(\hat{\mathbf{X}}|\mathbf{X})\sim\mathcal{N}(a_{i}\mathbf{x}_{0},b_{i}^{2}\mathbf{I})$,
the formula can approximate the denoised sample $\mathbf{x}_{0}$:

$\displaystyle\mathbb{E}[\mathbf{X}_{0}|{{\mathbf{X}_{i}}}]=({\mathbf{X}_{i}}+b_{i}^{2}\nabla_{{\mathbf{X}_{i}}}\log
p({\mathbf{X}_{i}}))/a_{i}.$ (2)

Here $a_{i}$ and $b_{i}$ are the drift and diffusion coefficients of the
forward diffusion process. Our approach uses this equation to approximate the
clean point cloud from the noisy samples at each time step during the reverse
diffusion process.

###  2.2 Point Cloud Smoothness Property

Point clouds have extensive applications in computer vision and robotics.
Nevertheless, raw point clouds are often affected by noise and irregularities
resulting from data acquisition techniques or environmental factors [23]. One
pivotal approach to improving the quality of point cloud generative models
involves implementing smoothness constraints, primarily utilizing graph
Laplacian matrices constructed from k-nearest neighbor (KNN) graphs [24]. This
method characterizes the point cloud’s smoothness by quantifying the average
difference between each point and its neighboring points, as expressed by the
equation:

$\displaystyle S(\mathbf{X})=\text{trace}(\mathbf{X}^{T}\mathbf{L}\mathbf{X})$
(3)

Here $\mathbf{X}$ represents the matrix of all points’ positions and
$\mathbf{L}$ stands for the Laplacian matrix derived from the KNN graph. Point
cloud generative models can attain smoother surfaces by incorporating this
smoothness constraint term into an optimization framework.

##  3 Approach

A point cloud consists of a set of points in 3D space. Each point cloud can be
denoted by
$\mathbf{X}=\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}\in\mathbb{R}^{N\times 3}$,
where $N$ is the number of points, and $\boldsymbol{x}_{i}$ is the 3D
coordinates for point $i$. The dataset comprises $M$ distinct point clouds,
denoted as $\mathcal{D}=\\{{\bf X}_{1},...{\bf X}_{M}\\}$. The goal is to
develop a model to generate a precise point cloud that faithfully captures the
datasets’ shape geometry and variability. We formulate point cloud generation
as a continuous particle diffusion process.

However, this problem introduces several unique challenges. Firstly, a single
object can have different shapes (e.g., various airplane categories).
Therefore, it is essential to model the global shape distribution and the
individual point cloud distribution. Secondly, the diffusion probabilistic
model captures the intermediate movement of points but may not consider local
geometric information within the point cloud, potentially leading to a non-
smooth surface during sampling. To address these challenges, we propose a
diffusion probabilistic model with a smoothness constraint to capture the
local geometric features within point clouds effectively. The overall model
structure is depicted in Figure1. The model consists of an encoder, a
diffusion decoder, and a latent diffusion module. The encoder
$q_{\mathbold{\phi}}(\bf z|\bf X)$ parameterized by $\mathbold{\phi}$ learns a
low dimensional latent embedding distribution for the given point cloud input
to encode its global shape feature. Given a shape embedding $\mathbf{z}$, the
diffusion decoder $p_{\boldsymbol{\psi}}(\bf X|\bf z)$ reconstructs the
original point cloud through a reverse diffusion process. Furthermore, to
accurately model the shape latent embedding $\bf z$ prior distribution, we use
a latent diffusion module $p_{{\mathbold{\theta}}}(\bf z)$ to transform a
standard normal distribution to the aggregated posterior of the encoder.

Figure 1: The overall structure of our model, containing three modules:
Encoder module $q_{\mathbold{\phi}}(\mathbf{z}|\mathbf{X})$; The latent
diffusion learns the latent prior distribution
$p_{\mathbold{\theta}}(\mathbf{z})$; The decoder diffusion module learns a
conditional diffusion model $p_{\boldsymbol{\psi}}{(\mathbf{X}|\mathbf{z})}$
to reconstruct the original shape $\mathbf{X}$ based on its embedding
$\mathbf{z}$. A smoothness constraint is incorporated during the decoder
reverse diffusion process.

###  3.1 Proposed Overall Model

Our model learns the point cloud distribution over a diverse shape pattern. As
indicated by [16], for the point cloud generation, it is not helpful to learn
the marginal distribution $p(\mathbf{X})$ across various shapes. Thereby, we
model the distribution into the factorized format and incorporate the latent
shape embedding distribution into consideration:
$p_{\mathbold{\theta},\boldsymbol{\psi}}(\mathbf{X},\mathbf{z})=p_{\boldsymbol{\psi}}(\mathbf{X}|\mathbf{z})p_{\mathbold{\theta}}(\mathbf{z})$,
where $\mathbf{z}\in\mathbb{R}^{d}$ is a latent d-dimensional variable that
embeds the global shape pattern of the point cloud,
$p_{\mathbold{\theta}}(\mathbf{z})$ represents the latent embedding prior
distribution. $p_{\boldsymbol{\psi}}(\mathbf{X}|\mathbf{z})$ is a conditional
diffusion model given a specific global shape condition $\mathbf{z}$. The
model trains the time-dependent score model
$s_{\boldsymbol{\psi}}(\mathbf{X},\mathbf{z},t)$ to learn the score function
of $p_{\boldsymbol{\psi}}(\mathbf{X}|\mathbf{z})$, which is the gradient of
log-likelihood $\nabla_{\mathbf{X}}\log
p_{\boldsymbol{\psi}}(\mathbf{X}|\mathbf{z})$. The global shape posterior is
approximated through another latent diffusion model
$q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{X})$. We maximize the Evidence Lower
Bound (ELBO) on the data log-likelihood as Equation 3.1 to optimize the
generative model.

$\displaystyle\max_{\boldsymbol{\phi},\boldsymbol{\psi},\mathbold{\theta}}~{}\mathcal{L}(\mathbold{\phi},\boldsymbol{\psi},\mathbold{\theta})$
$\displaystyle=\mathbb{E}_{q{\phi}(\mathbf{z}|\mathbf{X})}\big{[}\log{p_{\boldsymbol{\psi}}(\mathbf{X}|\mathbf{z})}\big{]}-\textrm{KL}\left(q{\phi}(\mathbf{z}|\mathbf{X})\;\|\;p{\theta}(\mathbf{z})\right)$
$\displaystyle=\mathbb{E}_{q{\phi}(\mathbf{z}|\mathbf{X})}\big{[}\log{p_{\boldsymbol{\psi}}(\mathbf{X}|\mathbf{z})}\big{]}+\mathbb{E}_{q{\phi}(\mathbf{z}|\mathbf{X})}\big{[}\log
p{\theta}(\mathbf{z})\big{]}+\mathbb{H}\big{[}q{\phi}(\mathbf{z}|\mathbf{X})\big{]}$
(4)

###  3.2 Sampling with Smoothness Constraint

The framework uses a conditioned diffusion probabilistic model to generate
point clouds with diverse shapes. However, it is still challenging to generate
the 3D point location distribution obeying the geometrical constraint on the
3D object.

We propose a graph Laplacian smoothness constraint on the decoder diffusion
model $p_{\boldsymbol{\psi}}({\bf X}|\mathbf{z})$. Specifically, we
incorporate the point cloud constraint into the sampling process. Suppose a
constraint $\mathbf{H}$ on the sample takes the form $p(\mathbf{H}|\bf{X})$,
we consider a probabilistic model from Bayesian perspective:

$\displaystyle
p(\bf{X}|\mathbf{H},\mathbf{z})=\frac{p(\mathbf{H}|\bf{X},\mathbf{z})p(\bf{X}|\mathbf{z})}{p(\mathbf{H},\mathbf{z})}$
(5)

We propose a Bayesian denoising method by taking the gradient of the
posterior:

$\displaystyle\nabla_{\bf{X}}\log
p(\bf{X}|\mathbf{H},\mathbf{z})=\nabla_{\bf{X}}\log
p(\bf{X}|\mathbf{z})+\nabla_{\bf{X}}\log p(\mathbf{H}|\bf{X},\mathbf{z})$ (6)

The first term is learned by the decoder diffusion score model
$s_{\boldsymbol{\psi}}(\mathbf{X},\mathbf{z},t)$, and the challenge is how to
add the second Bayesian constraints term to each time step of the reverse
diffusion process on noisy sample $\mathbf{X}_{t}$. To address this challenge,
We consider the constraints on the inferred clean point clouds based on
Tweedie’s formula:

$\displaystyle\hat{\mathbf{X}}=\frac{\mathbf{X}_{t}+b_{t}^{2}\mathbf{s}_{\boldsymbol{\psi}}(\mathbf{X},\mathbf{z},t)}{a_{t}}.$
(7)

Here $a_{t},b_{t}$ are the drift and diffusion coefficients of the forward
diffusion processes. Thus, we can construct a smoothness constraint based on
the estimated clean point cloud $\hat{\bf{X}}$. The constraint is proportional
to the graph Laplacian regularizer on the KNN graph of the point cloud to
encourage more uniform point distribution during the sampling process:

$\displaystyle\log
p(\mathbf{H}|{\mathbf{X}},\mathbf{z})\propto-\text{trace}(\hat{\mathbf{X}}^{T}\mathbf{L}\hat{\mathbf{X}}),\text{where}\
\hat{\mathbf{X}}={\mathbf{X}_{t}+b_{t}^{2}s_{\boldsymbol{\psi}}(\mathbf{X}_{t},\mathbf{z},t)}/{a_{t}}.$
(8)

Here $\hat{\bf{X}}$ is the estimated clean point cloud data at each reverse
diffusion time step, and $\bf{L}$ is the graph Laplacian matrix. Putting the
Bayesian denoising method together, the reverse diffusion under the smoothness
constraint can be represented by

$\displaystyle\mathbf{X}_{{t-1}}=\mathbf{X}_{{t}}+\underbrace{[f_{t}(\mathbf{X}_{t})-g_{t}^{2}\mathbf{s}_{\boldsymbol{\psi}}(\mathbf{X}_{t},\mathbf{z},t)]\Delta_{t}}_{\text{parameterized
by score model}}-\underbrace{\alpha\nabla_{\mathbf{X}_{t}}\\{\text{trace}\
(\hat{\mathbf{X}}_{t}^{T}\mathbf{L}\hat{\mathbf{X}}_{t})\\}}_{\text{smoothness
constraint}}+g_{t}(\bf{X}_{t})\mathbf{z},\mathbf{z}\sim\mathcal{N}(0,I).$ (9)

Here $\alpha$ controls the relative weight of the smoothness constraint, and
$f_{t},g_{t}$ are the drift and diffusion coefficients of the forward
diffusion process.

##  4 Experiment

This section describes the experimental results of the point cloud generation
task. The experiment aims to quantitatively compare our proposed framework
with several state-of-the-art models.

###  4.1 Experiment Setup

Baseline methods: We compare our proposed model against several SOTA models
for point cloud generation. Those models include r-GAN[19], l-GAN[18],
SoftFlow[25], and DPM[13].

Datasets: For our experiments, we utilized the ShapeNet dataset[26]. We follow
[12] to split dataset. Evaluation metrics: We adopt the following metrics:
Minimum Matching Distance(MMD), Coverage score(COV), 1-NN classifier accuracy
(1-NNA), and Relative Smoothness (RS). RS is computed by the graph Laplacian
as Eq.3. It is computed by:

$\displaystyle RS(P_{model},P_{data})=\|\mathbb{E}_{X\sim
p_{model}}\text{trace}(\mathbf{X}^{T}\mathbf{L}\mathbf{X})-\mathbb{E}_{\mathbf{X}\sim
p_{data}}\text{trace}(\mathbf{X}^{T}\mathbf{L}\mathbf{X})\|$ (10)

###  4.2 Generation Performance Comparison

Table 1: Comparison of shape generation on ShapeNet. MMD is multiplied by $10^{2}$. Numbers in bold indicate the best two methods Category | Model |  MMD ($10^{2},\downarrow$) |  COV (%,$\uparrow$) |  1-NNA (%,$\downarrow$)
---|---|---|---|---
Airplane | r-GAN | 2.31 | 14.32 | 96.79
l-GAN | 0.58 | 21.23 | 93.95
SoftFlow | 0.38 | 47.90 | 65.80
DPM | 0.57 | 33.83 | 86.91
Ours | 0.39 | 49.11 | 68.89
Ours+constraint | 0.38 | 49.63 | 66.75
Chair | r-GAN | 8.31 | 15.13 | 99.70
l-GAN | 2.01 | 29.31 | 99.70
SoftFlow | 1.68 | 47.43 | 60.05
DPM | 2.07 | 35.50 | 74.77
Ours | 1.61 | 47.89 | 61.15
Ours+constraint | 1.60 | 50.19 | 60.54
Car | r-GAN | 2.13 | 6.539 | 99.01
l-GAN | 1.23 | 23.58 | 88.78
SoftFlow | 0.86 | 44.60 | 60.09
DPM | 1.14 | 34.94 | 79.97
Ours | 0.85 | 47.16 | 71.88
Ours+constraint | 0.82 | 47.88 | 62.90

We compare our model with baseline models on the airplane, chair, and car
categories for point cloud generation. We consider two settings of our
proposed model. "Ours" represents the model without smoothness constraint and
the other with smoothness constraint. The evaluation results are summarized in
Table 1. The results show that our proposed model outperforms the GAN and
diffusion-based baseline models. Our model performs best in the MMD, COV, and
1-NNA measures. The model with constraint further improves the performance.
Thus, it demonstrated the effectiveness of the proposed constraint. We also
visualize generated samples in the Appendix.

Table 2: Comparison of smoothness performance on ShapeNet. RS is the relative smoothness Category | Model |  RS ($\downarrow$) | GT smoothness
---|---|---|---
Airplane | Our diffusion model | 64.90 | 1651.92
Our diffusion model + constraints | 37.52
Chair | Our diffusion model | 1041.35 | 4807.79
Our diffusion model + constraints | 983.92
Car | Our diffusion model | 203.02 | 2677.01
Our diffusion model + constraints | 177.90

###  4.3 Smoothness Performance Comparison

Furthermore, to demonstrate the performance of the constraint, we use the
relative smoothness (RS) metrics to evaluate. The results are shown in Table
2. The RS metrics measure the relative smoothness between sampled and ground
truth point clouds, and we also provide the ground-truth smoothness in the
last column. By incorporating the smoothness constraints into the generation
model, we can observe a reduction in the relative smoothness metrics. This
serves to validate the effectiveness of the proposed smoothness constraint in
generating a smoother surface.

##  5 Conclusions

In this paper, we introduce a novel diffusion-based model for point cloud
generation, incorporating not only the global geometry distribution but also
integrating a local smoothness constraint into the generation process.
Experimental results demonstrate that the proposed framework can produce
realistic and smoother point cloud samples. The limitation of our approach is
that we only consider the smoothness constraint on the point cloud.

## Acknowledgments and Disclosure of Funding

We thank anonymous reviewers for their valuable feedback. This work is
partially supported by the NSF CAREER Award #2239869.

## References

  * [1] Xiangyu Yue, Bichen Wu, Sanjit A Seshia, Kurt Keutzer, and Alberto L Sangiovanni-Vincentelli.  A lidar point cloud generator: from a virtual world to autonomous driving.  In Proceedings of the 2018 ACM on International Conference on Multimedia Retrieval, pages 458–464, 2018.
  * [2] Jiageng Mao, Xiaogang Wang, and Hongsheng Li.  Interpolated convolutional networks for 3d point cloud understanding.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 1578–1587, 2019.
  * [3] Qingyong Hu, Bo Yang, Linhai Xie, Stefano Rosa, Yulan Guo, Zhihua Wang, Niki Trigoni, and Andrew Markham.  Randla-net: Efficient semantic segmentation of large-scale point clouds.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 11108–11117, 2020.
  * [4] Yulan Guo, Hanyun Wang, Qingyong Hu, Hao Liu, Li Liu, and Mohammed Bennamoun.  Deep learning for 3d point clouds: A survey.  IEEE transactions on pattern analysis and machine intelligence, 43(12):4338–4364, 2020.
  * [5] François Pomerleau, Francis Colas, Roland Siegwart, et al.  A review of point cloud registration algorithms for mobile robotics.  Foundations and Trends® in Robotics, 4(1):1–104, 2015.
  * [6] Wenchong He, Zhe Jiang, Tingsong Xiao, Zelin Xu, Shigang Chen, Ronald Fick, Miles Medina, and Christine Angelini.  A hierarchical spatial transformer for massive point samples in continuous space.  arXiv preprint arXiv:2311.04434, 2023.
  * [7] Saining Xie, Jiatao Gu, Demi Guo, Charles R Qi, Leonidas Guibas, and Or Litany.  Pointcontrast: Unsupervised pre-training for 3d point cloud understanding.  In Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part III 16, pages 574–591. Springer, 2020.
  * [8] Wenchong He, Zhe Jiang, Chengming Zhang, and Arpan Man Sainju.  Curvanet: Geometric deep learning based on directional curvature for 3d shape analysis.  In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 2214–2224, 2020.
  * [9] Renrui Zhang, Ziyu Guo, Wei Zhang, Kunchang Li, Xupeng Miao, Bin Cui, Yu Qiao, Peng Gao, and Hongsheng Li.  Pointclip: Point cloud understanding by clip.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8552–8562, 2022.
  * [10] Hengshuang Zhao, Li Jiang, Chi-Wing Fu, and Jiaya Jia.  Pointweb: Enhancing local neighborhood features for point cloud processing.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 5565–5573, 2019.
  * [11] Qingyong Hu, Bo Yang, Guangchi Fang, Yulan Guo, Aleš Leonardis, Niki Trigoni, and Andrew Markham.  Sqn: Weakly-supervised semantic segmentation of large-scale 3d point clouds.  In European Conference on Computer Vision, pages 600–619. Springer, 2022.
  * [12] Guandao Yang, Xun Huang, Zekun Hao, Ming-Yu Liu, Serge Belongie, and Bharath Hariharan.  Pointflow: 3d point cloud generation with continuous normalizing flows.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 4541–4550, 2019.
  * [13] Shitong Luo and Wei Hu.  Diffusion probabilistic models for 3d point cloud generation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2837–2845, 2021.
  * [14] Ruojin Cai, Guandao Yang, Hadar Averbuch-Elor, Zekun Hao, Serge Belongie, Noah Snavely, and Bharath Hariharan.  Learning gradient fields for shape generation.  In European Conference on Computer Vision, pages 364–381. Springer, 2020.
  * [15] Yongbin Sun, Yue Wang, Ziwei Liu, Joshua Siegel, and Sanjay Sarma.  Pointgrow: Autoregressively learned point cloud generation with self-attention.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 61–70, 2020.
  * [16] Chun-Liang Li, Manzil Zaheer, Yang Zhang, Barnabas Poczos, and Ruslan Salakhutdinov.  Point cloud gan.  arXiv preprint arXiv:1810.05795, 2018.
  * [17] Dong Wook Shu, Sung Woo Park, and Junseok Kwon.  3d point cloud generative adversarial network based on tree structured graph convolutions.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 3859–3868, 2019.
  * [18] Diego Valsesia, Giulia Fracastoro, and Enrico Magli.  Learning localized generative models for 3d point clouds via graph convolution.  In International conference on learning representations, 2018.
  * [19] Panos Achlioptas, Olga Diamanti, Ioannis Mitliagkas, and Leonidas Guibas.  Learning representations and generative models for 3d point clouds.  In International conference on machine learning, pages 40–49. PMLR, 2018.
  * [20] Linqi Zhou, Yilun Du, and Jiajun Wu.  3d shape generation and completion through point-voxel diffusion.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5826–5835, 2021.
  * [21] Herbert E Robbins.  An empirical bayes approach to statistics.  In Breakthroughs in Statistics: Foundations and basic theory, pages 388–394. Springer, 1992.
  * [22] Kwanyoung Kim and Jong Chul Ye.  Noise2score: tweedie’s approach to self-supervised image denoising without clean images.  Advances in Neural Information Processing Systems, 34:864–874, 2021.
  * [23] Wenchong He and Zhe Jiang.  A survey on uncertainty quantification methods for deep neural networks: An uncertainty source perspective.  arXiv preprint arXiv:2302.13425, 2023.
  * [24] Jin Zeng, Gene Cheung, Michael Ng, Jiahao Pang, and Cheng Yang.  3d point cloud denoising using graph laplacian regularization of a low dimensional manifold model.  IEEE Transactions on Image Processing, 29:3474–3489, 2019.
  * [25] Hyeongju Kim, Hyeonseung Lee, Woo Hyun Kang, Joun Yeop Lee, and Nam Soo Kim.  Softflow: Probabilistic framework for normalizing flow on manifolds.  Advances in Neural Information Processing Systems, 33:16388–16397, 2020.
  * [26] Angel X Chang, Thomas Funkhouser, Leonidas Guibas, Pat Hanrahan, Qixing Huang, Zimo Li, Silvio Savarese, Manolis Savva, Shuran Song, Hao Su, et al.  Shapenet: An information-rich 3d model repository.  arXiv preprint arXiv:1512.03012, 2015.
  * [27] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  arXiv preprint arXiv:2011.13456, 2020.
  * [28] Arash Vahdat, Karsten Kreis, and Jan Kautz.  Score-based generative modeling in latent space.  Advances in Neural Information Processing Systems, 34:11287–11302, 2021.
  * [29] Xiaohui Chen, Yukun Li, Aonan Zhang, and Li-ping Liu.  Nvdiff: Graph generation through the diffusion of node vectors.  arXiv preprint arXiv:2211.10794, 2022.
  * [30] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 652–660, 2017.
  * [31] Lars Mescheder, Michael Oechsle, Michael Niemeyer, Sebastian Nowozin, and Andreas Geiger.  Occupancy networks: Learning 3d reconstruction in function space.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 4460–4470, 2019.

## Appendix

##  6 Background: Diffusion-based Generative Model

We consider the diffusion process over data $\mathbf{X}$ indexed with time
$t\in[0,1],\\{\mathbf{X}_{t}\\}_{t=0}^{1}$. Here $\mathbf{X}_{0}\sim
p_{0}(\mathbf{X})=p_{data}$, $\mathbf{X}_{1}\sim p_{1}(\mathbf{X})$ , where
$p_{1}(\mathbf{X})$ approximates the Gaussian distribution
$\mathcal{N}(\mathbf{0},\mathbf{I})$. The diffusion process can be modeled by
Itô SDE[27]:

$\displaystyle d\mathbf{X}_{t}=f_{t}(\mathbf{X},t)dt+g_{t}d\mathbf{w},$ (11)

$f_{t}$ is a drift coefficient, and $g_{t}$ is a diffusion coefficient. The
forward diffusion process smoothly transforms $\mathbf{X}$ by adding
infinitesimal noise $dw$ at each infinitesimal step $dt$.

To generate samples following the data distribution, we start from the prior
distribution and follow the reverse diffusion process, which is defined by:

$\displaystyle
d\mathbf{X}_{t}=[f_{t}(\mathbf{X}_{t})-g_{t}^{2}\nabla_{\mathbf{X}_{t}}\log
p_{t}(\mathbf{X}_{t})]dt+g_{t}d\bar{\mathbf{w}},$ (12)

The diffusion model training process learns the time-dependent score model
$\mathbf{s}_{\theta}(\mathbf{X},t)$ and minimizes the least square loss
between the score function and the model output

$\displaystyle\min\int_{0}^{1}\mathbb{E}_{p_{t}(\mathbf{X})}[\|\mathbf{s}_{\theta}(\mathbf{X},t)-\nabla_{\mathbf{X}}\log
p_{t}(\mathbf{X})\|^{2}]dt$ (13)

Thus, the sampling process replace the unknown score function
$\nabla_{\mathbf{X}_{t}}\log p_{t}(\mathbf{X}_{t})$ in Eq 12 with the score
model $\mathbf{s}_{\theta}(\mathbf{X},t)$.

##  7 Details on Approach

###  7.1 Diffusion-based Prior Model

The prior model aims to capture the distribution of the global geometry of all
point cloud shapes. A common choice is to use a fixed prior (e.g.,
$\mathbf{z}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$) in variational auto-
encoder, but it can produce holes, then there exists a mismatch between the
aggregated posterior
$\sum_{i}q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{X}^{(i)})$ and the prior
$p(\mathbf{z})$. Then, this can lead to some regions in the latent space that
cannot be decoded into data-like samples.

Our proposed model uses a learnable prior model $p_{\theta}(\mathbf{z})$. The
prior is flexible, such that it can catch up with the aggregated posterior,
and then they can match each other better, thus alleviating the problem of
prior holes. Here, we propose to use the diffusion process to learn the prior
over the shape distribution. Specifically, given a global geometry encoding
$\mathbf{z}_{0}\sim q_{\boldsymbol{\phi}}(\mathbf{z}_{0}|\mathbf{X})$, we
diffused $z_{0}$ through a SDE defined in Eq. 11 to the standard Normal
distribution $p(\mathbf{z}_{1})=\mathcal{N}(\mathbf{0},\mathbf{I})$, and we
train a time-dependent score model
$\mathbf{s}_{\mathbold{\theta}}(\mathbf{z}_{t},t)$ to approximate the true
score $\nabla_{\mathbf{z}_{t}}\log p_{\mathbold{\theta}}(\mathbf{z}_{t})$. For
the generation process, the reversed SDE based on Eq. 12 is employed to
generate samples that follow the prior distribution
$p_{\mathbold{\theta}}(\mathbf{z}_{0})$. We use several layers of 1D
convolution with residual connection for the score model. The details are in
the experiment section.

###  7.2 Diffusion-based Conditional Decoder

In this section, we introduce the framework of the conditional diffusion model
in the decoder part. Here, we consider the point cloud generation as a
conditional diffusion process, where the condition is the global geometry
encoding $\mathbf{z}\sim q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{X})$.
Similarly, we diffuse the points of the point cloud (particles) to the
standard Normal distribution and learn a conditional score model
$\mathbf{s}_{\boldsymbol{\psi}}(\mathbf{X}_{t},\mathbf{z},t)$ to approximate
the reverse moving direction of the points during the forward diffusion
process (equivalent to the score $\nabla_{\mathbf{X}_{t}}\log
p_{\boldsymbol{\psi}}(\mathbf{X}_{t}|\mathbf{z})$). For sampling, given a
global geometry condition $\mathbf{z}$, we can use the score model
$s_{\boldsymbol{\psi}}(\mathbf{X}_{t},\mathbf{z},t)$ through the reverse SDE
to move the noisy point cloud to the original point cloud distribution.

Specifically, the input of this network is the perturbed point cloud with
Gaussian noise sampled from $\mathbf{X}_{t}\sim
p_{t}(\mathbf{X}_{t}|\mathbf{X})$ and the time embedding. We also concatenate
the latent code $\mathbf{z}$ as the condition to the input of each layer.

###  7.3 Training algorithm

In the training objective of Eq.3.1, we estimate the entropy term
$\mathbb{H}\big{[}q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{X})\big{]}$ via
Monte Carlo samples. And our decoder is a continuous diffusion-based
conditional decoder, then the reconstruction term
$\mathbb{E}_{q_{\boldsymbol{\phi}}}(\mathbf{z}|\mathbf{X})\big{[}\log{p_{\boldsymbol{\psi}}(\mathbf{X}|\mathbf{z})}\big{]}$
is approximated by mean square loss of the score function

$\displaystyle\mathbb{E}_{q_{\boldsymbol{\phi}}}(\mathbf{z}|\mathbf{X})\big{[}\log{p_{\boldsymbol{\psi}}(\mathbf{X}|\mathbf{z})}\big{]}$
$\displaystyle=\mathbb{E}_{t,p_{t}(\mathbf{X})}\big{[}\frac{g(t)^{2}}{2}\|\mathbf{s}_{\boldsymbol{\psi}}(\mathbf{X}_{t},\mathbf{z},t)-\nabla_{\mathbf{X}}\log
p_{t}(\mathbf{X}_{t}|\mathbf{z})\|^{2}_{2}\big{]}\big{]}$ (14)

and we optimize the cross-entropy term through score matching following [28,
29] with $t\sim U[0,1]$:

$\displaystyle\mathbb{E}_{q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{X})}$
$\displaystyle\big{[}\log
p_{\theta}(\mathbf{z})\big{]}=\mathbb{E}_{t,q_{\mathbold{\phi}}(\mathbf{z}_{0},\mathbf{z}_{t}|\mathbf{X})}\Big{[}\frac{g(t)^{2}}{2}\|\mathbf{s}{\theta}(\mathbf{z}_{t},t)-\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}_{t})\|^{2}_{2}\Big{]}$ (15)

The training process is described in Algorithm 1

Algorithm 1 Training Procedure

Input: point cloud training data $p_{data}(\mathbf{X})$

for epoch = $0$ to $n$ do

Sample $\mathbf{X}_{0}\sim p_{data}(X)$, $t\sim\text{Uniform}(0,1)$

Encoder:
$\mathbf{z}_{0}=q_{\mathbold{\phi}}(\mathbf{z}|\mathbf{X}_{0}),\mathbf{z}_{t}\sim
p_{t}(\mathbf{z}_{t}|\mathbf{z}_{0})$

Latent diffusion loss:
$L_{\mathbf{z}}^{t}\leftarrow\frac{g(t)^{2}}{2}\|\mathbf{s}{\theta}(\mathbf{z}_{t},t)-\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}_{t})\|^{2}_{2}$

Sample $\mathbf{X}_{t}\sim p_{t}(\mathbf{X}_{t}|\mathbf{X}_{0})$

Recon loss:
$L_{\mathbf{x}}^{t}\leftarrow\frac{g(t)^{2}}{2}\|\mathbf{s}_{\boldsymbol{\psi}}(\mathbf{X}_{t},\mathbf{z},t)-\nabla_{\mathbf{X}}\log
p_{t}(\mathbf{X}_{t}|\mathbf{z})\|^{2}_{2}$

Apply stochastic gradient on the neural network parameters with
$\nabla_{\mathbold{\theta}}(L_{\mathbf{z}}^{t}+L_{\mathbf{x}}^{t}+\mathbb{H}[q{\phi}(\mathbf{z}|\mathbf{X}_{t})])$

end for

Algorithm 2 Sampling Procedure

Input: Latent score model $\mathbf{s}{\theta}(\mathbf{z}_{t},t)$ and decoder
score model $\mathbf{s}{\phi}(\mathbf{X}_{t},\mathbf{z},t)$

Sample
$\mathbold{\epsilon},\mathbf{x}_{\text{init}}\sim\mathcal{N}(\mathbold{\epsilon};\mathbf{0},\mathbf{I})$

Latent sampling: $\mathbold{\epsilon}\xrightarrow{\text{diffusion}}\mathbf{z}$

$\mathbf{X}_{t_{n}}=\mathbf{X}_{\text{init}}$, $t=n$, $\Delta_{t}=\frac{1}{n}$

for $i=n$ to $1$ do

$\small\mathbf{X}_{{t-1}}=\mathbf{X}_{{t}}+{[f_{t}(\mathbf{X}_{t})-g_{t}^{2}\mathbf{s}_{\boldsymbol{\psi}}(\mathbf{X}_{t},\mathbf{z},t)]\Delta_{t}}-{\alpha\nabla_{\mathbf{X}_{t}}\\{\text{trace}\
(\hat{\mathbf{X}}_{t}^{T}\mathbf{L}\hat{\mathbf{X}}_{t})\\}}+g_{t}(\bf{X}_{t})\mathbf{z},\
\mathbf{z}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$

end for

return $\mathbf{X}_{0}$

##  8 More Details on Experiments

###  8.1 Experiment Setup

Datasets: In our experimental evaluation, we employed the widely used ShapeNet
dataset[26], a comprehensive benchmark for 3D shape understanding. The dataset
comprises 51,127 shapes distributed across 55 different categories. To
maintain consistency with prior work and ensure fair comparisons, we adopted
the same dataset split as presented in [12].

Within the ShapeNet dataset, e each 3D shape is depicted as a point cloud
consisting of 2048 points uniformly sampled on the object’s surfaces. This
representation of point clouds effectively captures the intricate geometric
details of the shapes, making it a fitting input for our 3D shape
understanding tasks.

For our comparative analysis, we specifically focused on three categories:
airplanes, chairs, and cars. These categories represent different types of
objects with varying shapes and complexities, making them suitable for
assessing the performance of our model under diverse scenarios.

###  8.2 Implementation Detail

Hyperparameters: The latent code dimension is set to 256. For the training
stage, we set the batch size as 32. We utilized Adam optimizer with learning
rate $2\times 10^{-3}$ for the encoder and learning rate $2\times 10^{-4}$ for
the diffusion decoder module, $1\times 10^{-4}$ for the latent diffusion
module. The learning rate was kept constant for the first 1000 epochs and
decreased linearly from 1000 to 2000 epochs. For sampling, We sample each
point cloud with N = 2048 points. The time step is set to be $1000$ steps for
sampling.

Model architecture: We use PointNet [30] as an encoder to map the input point
cloud into a 512-dimension latent feature code $\mathbf{z}$. For the decoder
and the prior diffusion module, we utilize OccNet [31]. We stacked 6 ResNet
blocks with 256 dimensions for every hidden layer.

Evaluation metrics:

Suppose $P_{r}$ is the set of reference point clouds for the generation task,
and $P_{g}$ is the set of generated point clouds. We adopt the following
metrics:

  1. 1.

Minimum Matching Distance(MMD)

$\footnotesize MMD(P_{r},P_{g})=\frac{1}{|P_{r}|}\sum_{P\in P_{r}}\min_{Q\in
P_{g}}D(P,Q)$

MMD calculates the average of the closest distances to each reference point
cloud to measure the quality of the generated samples.

  2. 2.

Coverage score(COV)

$\displaystyle\footnotesize COV(P_{r},P_{g})=\frac{|\\{\arg min_{P\in
P_{r}}D(P,Q)|Q\in P_{g}\\}|}{|P_{r}|}$

COV measures the variety of the generated samples to detect mode collapse.

  3. 3.

1-NN classifier accuracy (1-NNA): The score is computed using a 1-NN
classifier to classify the generated samples and the reference point cloud.
When the quality of the generated samples is high, the 1-NNA score should be
close to 50%.

  4. 4.

Relative Smoothness (RS): The smoothness of the sampled point cloud is
computed by the graph Laplacian as Eq.3. It is computed by:

$\displaystyle RS(P_{model},P_{data})=\|\mathbb{E}_{X\sim
p_{model}}\text{trace}(\mathbf{X}^{T}\mathbf{L}\mathbf{X})-\mathbb{E}_{\mathbf{X}\sim
p_{data}}\text{trace}(\mathbf{X}^{T}\mathbf{L}\mathbf{X})\|$

Figure 2: Visualization of some generated samples by our model Figure 3:
Sensitivity analysis on varying $K$ in KNN graph construction for airplane
dataset

In addition, we perform a sensitivity analysis to investigate how the choice
of neighbors of the KNN graph influences the smoothness performance.
Specifically, we conduct experiments by varying $K$ from $5$ to $35$ to
construct the KNN graph on the airplane datasets. We observe that the sample
smoothness is constantly lower than that generated by the diffusion model
without smoothness constraint ($1715$), and the best result is achieved when
$K=30$.


IEEEexample:BSTcontrol

# Speech Enhancement by Noise Self-Supervised Rank-Constrained Spatial
Covariance Matrix Estimation via Independent Deeply Learned Matrix Analysis

Sota Misawa1, Norihiro Takamune1, Tomohiko Nakamura1, Daichi Kitamura2,
Hiroshi Saruwatari1,
Masakazu Une3 and Shoji Makino34 1 The University of Tokyo, Tokyo, Japan
2 National Institute of Technology, Kagawa College, Kagawa, Japan 3 University
of Tsukuba, Ibaraki, Japan 4 Waseda University, Fukuoka, Japan

###### Abstract

Rank-constrained spatial covariance matrix estimation (RCSCME) is a method for
the situation that the directional target speech and the diffuse noise are
mixed. In conventional RCSCME, independent low-rank matrix analysis (ILRMA) is
used as the preprocessing method. We propose RCSCME using independent deeply
learned matrix analysis (IDLMA), which is a supervised extension of ILRMA. In
this method, IDLMA requires deep neural networks (DNNs) to separate the target
speech and the noise. We use Denoiser, which is a single-channel speech
enhancement DNN, in IDLMA to estimate not only the target speech but also the
noise. We also propose noise self-supervised RCSCME, in which we estimate the
noise-only time intervals using the output of Denoiser and design the prior
distribution of the noise spatial covariance matrix for RCSCME. We confirm
that the proposed methods outperform the conventional methods under several
noise conditions.

##  1 Introduction

In this study, we deal with speech enhancement. In particular, we deal with
the problem of extracting directional target speech from diffuse noise. This
situation arises, for example, when the speaker is close to the microphones in
a noisy environment such as a crowded place or a train station.

From a practical perspective, it is effective to use blind source separation
(BSS) [1], which does not require any spatial information or characteristics
of each source. In a determined or overdetermined situation (number of
microphones $M$ $\geq$ number of sources $N$), independent vector analysis
(IVA) [2, 3, 4] and independent low-rank matrix analysis (ILRMA) [5] have been
proposed for BSS. ILRMA assumes that the power spectrogram for each source is
modeled by nonnegative matrix factorization (NMF) [6]. Since these methods
assume that the rank of the spatial covariance matrix (SCM) of each source is
unity, we call these methods rank-1 methods. In a rank-1 method, the noise
remains in the separated target signal when the noise has a full-rank SCM,
such as diffuse noise [7, 8, 9].

Multichannel NMF (MNMF) [10, 11] and FastMNMF [12, 13] have been proposed as
methods to model full-rank SCMs for each source. MNMF adopts the full-rank
nature for the source model [14] and can handle diffuse sources, but it is
computationally expensive and highly dependent on initial values [5]. In
FastMNMF, the computational complexity is reduced by assuming that the SCMs
are jointly diagonalizable. However, the problem that the performance is not
robust to changes in the initial value still remains.

Rank-constrained spatial covariance matrix estimation (RCSCME) [15, 16] has
been proposed as a method of blind speech enhancement for mixed signals of
directional target speech and diffuse noise. RCSCME uses a rank-1 method such
as ILRMA in advance. Since a rank-1 method accurately obtains $M-1$ noise
signals, we can obtain the rank-($M-1$) component of the noise SCM from these
signals. RCSCME estimates the time-varying variance of the target speech and
the diffuse noise, and complements the deficient rank-1 component of the noise
SCM at the same time. As shown in [15], this method is robust and converges
rapidly because there are fewer estimated parameters than in MNMF and
FastMNMF.

For the case that training data are available, many supervised source
separation methods using a deep neural network (DNN) have been proposed. For a
determined or overdetermined situation, independent deeply learned matrix
analysis (IDLMA) [17] has been proposed to improve source separation
performance by changing the source model of ILRMA from NMF to DNNs. Owing to
the appropriately trained DNNs, IDLMA achieves higher performance than ILRMA.
Since IDLMA is a rank-1 method similarly to ILRMA, the noise remains in the
separated target signal in a situation with diffuse noise.

In this paper, we propose the use of IDLMA, which is a rank-1 method, as the
preprocessing method of RCSCME. Since IDLMA was originally proposed for
musical source separation, we should prepare DNNs that estimate the power
spectrogram of speech or noise. One of the DNN-based single-channel speech
enhancement methods is Denoiser [18], which estimates the target speech from a
mixed-signal input in the waveform domain. Whereas IDLMA requires DNNs for
each source, Denoiser extracts only the target speech. Thus, we also propose a
scheme to estimate the power spectrogram of the noise using Denoiser.
Furthermore, considering the high capability of Denoiser to reduce noise, we
can find the noise-only time intervals of the observed signal using the output
of Denoiser. The SCM of this noise-only signal is a good approximation of the
noise SCM. Thus, we design the prior distribution of the noise using this
approximated noise SCM for RCSCME, and we call this method noise self-
supervised RCSCME. We derive the new update rule on the basis of the EM
algorithm for noise self-supervised RCSCME. We conduct a simulated experiment
to verify the effectiveness of the proposed methods.

##  2 Conventional Methods

###  2.1 ILRMA and IDLMA

Let $\bm{s}_{ij}=(s_{ij,1},\ldots,s_{ij,N})^{\mathsf{T}}$,
$\bm{x}_{ij}=(x_{ij,1},\ldots,x_{ij,M})^{\mathsf{T}}$, and
$\bm{y}_{ij}=(y_{ij,1},\ldots,y_{ij,N})^{\mathsf{T}}$ denote the short-time
Fourier transform (STFT) of the source, observed, and separated signals,
respectively, where $i=1,\ldots,I$ and $j=1,\ldots,J$ are the indices of the
frequency bin and time frame, respectively. Here, T denotes transpose. We also
represent the spectrograms for each signal as
$S_{n},X_{m},Y_{n}\in\mathbb{C}^{I\times J}$ whose $(i,j)$th elements are
$s_{ij,n},x_{ij,m}$, and $y_{ij,n}$, where $n=1,\ldots,N$ and $m=1,\ldots,M$
are the source and channel indices, respectively.

When the reverberation time is sufficiently shorter than the window length of
the STFT and each source is a point source, we can denote the observed signal
as

$\displaystyle\bm{x}_{ij}=A_{i}\bm{s}_{ij},$ (1)

where $A_{i}=\quantity(\bm{a}_{i,1},\ldots,\bm{a}_{i,N})\in\mathbb{C}^{M\times
N}$ is a mixing matrix. Assuming that $N=M$ and $A_{i}$ is invertible, we can
obtain the separated signal by estimating its inverse matrix
$W_{i}=\quantity(\bm{w}_{i,1},\ldots,\bm{w}_{i,N})^{\mathsf{H}}\in\mathbb{C}^{N\times
M}$, which is called the demixing matrix, such that

$\displaystyle\bm{y}_{ij}=W_{i}\bm{x}_{ij},$ (2)

where H denotes the Hermitian transpose.

In ILRMA[5] and IDLMA[17], it is assumed that the separated signals are
mutually independent and that the $n$th separated signal in each time-
frequency frame is generated by the univariate complex Gaussian distribution
as

$\displaystyle
y_{ij,n}\sim\mathcal{N}_{\mathrm{c}}\left(0,\sigma_{ij,n}^{2}\right),$ (3)

where $0$ and $\sigma_{ij,n}^{2}$ are the mean and variance, respectively. We
use the negative log-likelihood of the observed signal as the cost function:

$\displaystyle\mathcal{L}_{\mathrm{I}}=$
$\displaystyle-\log\prod_{i,j,n}p(x_{ij,n})$ $\displaystyle=$
$\displaystyle-\log\prod_{i,j,n}p(y_{ij,n})-2J\sum_{i}\log|\det W_{i}|$
$\displaystyle=$
$\displaystyle\sum_{i,j,n}\quantity(\frac{|w_{in}^{\mathsf{H}}x_{ij}|^{2}}{\sigma_{ij,n}^{2}}+2\log\sigma_{ij,n})$
$\displaystyle-2J\sum_{i}\log|\det W_{i}|+\mathrm{const.},$ (4)

where $\mathrm{const.}$ includes the terms independent of the target
variables. In ILRMA, $\sigma_{ij,n}^{2}$ is modeled by NMF [6] as

$\displaystyle\sigma_{ij,n}^{2}=\sum_{k}t_{ik,n}v_{kj,n},$ (5)

where $t_{ik,n}\geq 0$ and $v_{kj,n}\geq 0$ are the NMF variables. Instead, in
IDLMA, we use a source separation DNN that estimates the power spectrogram:

$\displaystyle\sigma_{ij,n}^{2}$
$\displaystyle=\max\quantity{\absolutevalue{\quantity(\mathrm{DNN}_{n}\quantity(Y_{n}))_{ij}}^{2},\varepsilon},$
(6)

where $\varepsilon>0$ is a small number used to avoid numerical instability.
Here, $\mathrm{DNN}_{n}$ is the trained DNN corresponding to the $n$th source
and its input is the current separated signal.

Since diffuse noise is not a point source, ILRMA and IDLMA cannot separate
noise in the same direction as the directional target speech in principle [7,
8, 9]. Therefore, we cannot directly apply ILRMA and IDLMA to diffuse noise
reduction, which often arises in practical speech enhancement application.

###  2.2 RCSCME

RCSCME [15, 16] deals with the problem of speech enhancement when directional
target speech and diffuse noise are mixed. In [15], it is assumed that the
generative model of the observed signal is the multivariate complex Gaussian
distribution

$\displaystyle
p(x_{ij};0,R_{ij}^{(\mathrm{o})})\propto\frac{\exp(-x_{ij}^{\mathsf{H}}R_{ij}^{(\mathrm{o})}x_{ij})}{\det
R_{ij}^{(\mathrm{o})}},$ (7)

where $R_{ij}^{(\mathrm{o})}$ denotes the covariance matrix of the observed
signal. $R_{ij}^{(\mathrm{o})}$ is modeled as

$\displaystyle
R_{ij}^{(\mathrm{o})}=r_{ij}^{(\mathrm{t})}R_{i}^{(\mathrm{t})}+r_{ij}^{(\mathrm{n})}R_{i}^{(\mathrm{n})},$
(8)

where $r_{ij}^{(\mathrm{t})}>0$ and $r_{ij}^{(\mathrm{n})}>0$ are the
variances of the target speech and the noise, and $R_{i}^{(\mathrm{t})}$ and
$R_{i}^{(\mathrm{n})}$ are the SCMs of the target speech and the noise,
respectively. To induce the sparsity of the target speech, we assume the
inverse gamma distribution for the prior distribution on
$r_{ij}^{(\mathrm{t})}$,

$\displaystyle
p(r_{ij}^{(\mathrm{t})};\alpha,\beta)\propto\quantity(r_{ij}^{(\mathrm{t})})^{-\alpha-1}\exp(-\frac{\beta}{r_{ij}^{(\mathrm{t})}}),$
(9)

where $\alpha>0$ and $\beta>0$ are the shape parameter and scale parameter,
respectively. Using the estimated mixing matrix $W_{i}^{-1}=A_{i}$, we model
the SCMs of the target source and noise as

$\displaystyle R_{i}^{(\mathrm{t})}$
$\displaystyle=\bm{a}_{i,n_{t}}\bm{a}_{i,n_{t}}^{\mathsf{H}},$ (10)
$\displaystyle R_{i}^{(\mathrm{n})}$
$\displaystyle={{R^{\prime}}_{i}^{(\mathrm{n})}}+\lambda_{i}\bm{v}_{i}\bm{v}_{i}^{\mathsf{H}},$
(11) $\displaystyle{{R^{\prime}}_{i}^{(\mathrm{n})}}$
$\displaystyle=\frac{1}{J}\sum_{j}\hat{y}_{ij}^{(\mathrm{n})}\quantity(\hat{y}_{ij}^{(\mathrm{n})})^{\mathsf{H}},$
(12) $\displaystyle\hat{y}_{ij}^{(\mathrm{n})}$
$\displaystyle=A_{i}\left(w_{i,1}^{\mathsf{H}}x_{ij},\ldots,w_{i,n_{t}-1}^{\mathsf{H}}x_{ij},0,\right.$
$\displaystyle\qquad\qquad\quad\left.w_{i,n_{t}+1}^{\mathsf{H}}x_{ij},\ldots,w_{i,N}x_{ij}\right)^{\mathsf{T}},$
(13)

where $n_{t}$ denotes the target index and
${{R^{\prime}}_{i}^{(\mathrm{n})}}\in\mathbb{C}^{M\times M}$ is the
rank-($M-1$) SCM of the diffuse noise estimated by a rank-1 method. Since the
rank of ${{R^{\prime}}_{i}^{(\mathrm{n})}}$ is $M-1$, $\bm{v}_{i}$ must be
linearly independent of any column vectors of
${{R^{\prime}}_{i}^{(\mathrm{n})}}$ in order that $R_{i}^{(\mathrm{n})}$ has
full rank. For example, the unit eigenvector corresponding to the zero
eigenvalue of ${{R^{\prime}}_{i}^{(\mathrm{n})}}$ satisfies this condition for
$\bm{v}_{i}$. RCSCME simultaneously estimates $\lambda_{i}$,
$r_{ij}^{(\mathrm{t})}$, and $r_{ij}^{(\mathrm{n})}$ by maximum a posteriori
estimation. Using (7) and (9), we calculate the log-posterior of the observed
signal as

$\displaystyle\mathcal{L}_{\mathrm{RC}}(\Theta_{\mathrm{c}})=$
$\displaystyle\log\prod_{i,j}p(\bm{x}_{ij}|\Theta_{\mathrm{c}})p(r_{ij}^{(\mathrm{t})};\alpha,\beta)$
$\displaystyle=$
$\displaystyle-\sum_{i,j}\Biggl{(}\bm{x}_{ij}^{\mathsf{H}}R_{ij}^{(\mathrm{o})}\bm{x}_{ij}+\log\det
R_{ij}^{(\mathrm{o})}\Biggr{.}$ $\displaystyle+\left.(\alpha+1)\log
r_{ij}^{(\mathrm{t})}+\frac{\beta}{r_{ij}^{(\mathrm{t})}}\right)+\mathrm{const.},$
(14)

where
$\Theta_{\mathrm{c}}=\\{\lambda_{i},r_{ij}^{(\mathrm{t})},r_{ij}^{(\mathrm{n})}\\}$.
We maximize the log-posterior function using the EM algorithm with the
following $Q$ function (see [15] for details):

$\displaystyle
Q_{\mathrm{RC}}\quantity(\Theta_{\mathrm{c}};\tilde{\Theta}_{\mathrm{c}})=$
$\displaystyle-\sum_{i,j}\Bigg{(}(\alpha+2)\log
r_{ij}^{(\mathrm{t})}+\frac{\hat{r}_{ij}^{(\mathrm{t})}+\beta}{r_{ij}^{(\mathrm{t})}}$
$\displaystyle+M\log r_{ij}^{(\mathrm{n})}+\log\det R_{i}^{(\mathrm{n})}$
$\displaystyle+\frac{\tr(\hat{R}_{ij}^{(\mathrm{n})}\quantity(R_{i}^{(\mathrm{n})})^{-1})}{r_{ij}^{(\mathrm{n})}}\Bigg{)}+\mathrm{const.},$
(15)

where
$\tilde{\Theta}_{\mathrm{c}}=\\{\tilde{\lambda}_{i},\tilde{r}_{ij}^{(\mathrm{t})},\tilde{r}_{ij}^{(\mathrm{n})}\\}$
is the set of the up-to-date parameters, and $\hat{r}_{ij}^{(\mathrm{t})}$ and
$\hat{R}_{ij}^{(\mathrm{n})}$ are calculated in the E-step as follows:

$\displaystyle\tilde{R}_{i}^{(\mathrm{n})}=$
$\displaystyle{{R^{\prime}}_{i}^{(\mathrm{n})}}+\tilde{\lambda}_{i}\bm{v}_{i}\bm{v}_{i}^{\mathsf{H}},$
(16) $\displaystyle R_{ij}^{(\mathrm{o})}=$
$\displaystyle\tilde{r}_{ij}^{(\mathrm{t})}R_{i}^{(\mathrm{t})}+\tilde{r}_{ij}^{(\mathrm{n})}\tilde{R}_{i}^{(\mathrm{n})},$
(17) $\displaystyle\hat{r}_{ij}^{(\mathrm{t})}=$
$\displaystyle\tilde{r}_{ij}^{(\mathrm{t})}-\quantity(\tilde{r}_{ij}^{(\mathrm{t})})^{2}\bm{a}_{i,n_{t}}^{\mathsf{H}}\quantity(R_{ij}^{(\mathrm{o})})^{-1}\bm{a}_{i,n_{t}}$
$\displaystyle+\quantity|\tilde{r}_{ij}^{(\mathrm{t})}\bm{x}_{ij}^{\mathsf{H}}\quantity(R_{ij}^{(\mathrm{o})})^{-1}\bm{a}_{i,n_{t}}|^{2},$
(18) $\displaystyle\hat{R}_{ij}^{(\mathrm{n})}=$
$\displaystyle\tilde{r}_{ij}^{(\mathrm{n})}\tilde{R}_{i}^{(\mathrm{n})}-\quantity(\tilde{r}_{ij}^{(\mathrm{t})})^{2}\tilde{R}_{i}^{(\mathrm{n})}\quantity(R_{ij}^{(\mathrm{o})})^{-1}\tilde{R}_{i}^{(\mathrm{n})}$
$\displaystyle+\quantity(\tilde{r}_{ij}^{(\mathrm{t})})^{2}\tilde{R}_{i}^{(\mathrm{n})}\quantity(R_{ij}^{(\mathrm{o})})^{-1}\bm{x}_{ij}\bm{x}_{ij}^{\mathsf{H}}\quantity(R_{ij}^{(\mathrm{o})})^{-1}\tilde{R}_{i}^{(\mathrm{n})}.$
(19)

In the M-step, by finding the stationary point of $Q_{\mathrm{RC}}$, we obtain
the update rules as

$\displaystyle r_{ij}^{(\mathrm{t})}$
$\displaystyle\leftarrow\frac{\hat{r}_{ij}^{(\mathrm{t})}+\beta}{\alpha+2},$
(20) $\displaystyle\lambda_{i}$
$\displaystyle\leftarrow\bm{u}_{i}^{\mathsf{H}}\quantity(\frac{1}{J}\sum_{j}\frac{1}{\tilde{r}_{ij}^{(\mathrm{n})}}\hat{R}_{ij}^{(\mathrm{n})})\bm{u}_{i},$
(21) $\displaystyle R_{i}^{(\mathrm{n})}$
$\displaystyle\leftarrow{{R^{\prime}}_{i}^{(\mathrm{n})}}+{\lambda}_{i}\bm{v}_{i}\bm{v}_{i}^{\mathsf{H}},$
(22) $\displaystyle r_{ij}^{(\mathrm{n})}$
$\displaystyle\leftarrow\frac{1}{M}\tr(\hat{R}_{ij}^{(\mathrm{n})}\quantity(R_{i}^{(\mathrm{n})})^{-1}),$
(23)

where $\bm{u}_{i}$ is the eigenvector corresponding to the zero eigenvalues of
${{R^{\prime}}_{i}^{(\mathrm{n})}}$ satisfying
$\bm{u}_{i}^{\mathsf{H}}\bm{v}_{i}=1$. Here, we use the following property for
$R_{i}^{(\mathrm{n})}={{R^{\prime}}_{i}^{(\mathrm{n})}}+\lambda_{i}\bm{v}_{i}\bm{v}_{i}^{\mathsf{H}}$,
as proved in [16], in obtaining the update rule of $\lambda_{i}$:

$\displaystyle\log\det R_{i}^{(\mathrm{n})}$
$\displaystyle=\log\lambda_{i}+\mathrm{const.},$ (24)
$\displaystyle\quantity(R_{i}^{(\mathrm{n})})^{-1}$
$\displaystyle=\quantity(E-\bm{u}_{i}\bm{v}_{i}^{\mathsf{H}})\quantity({{R^{\prime}}_{i}^{(\mathrm{n})}})^{+}\quantity(E-\bm{v}_{i}\bm{u}_{i}^{\mathsf{H}})+\frac{1}{\lambda_{i}}\bm{u}_{i}\bm{u}_{i}^{\mathsf{H}},$
(25)

where $E\in\mathbb{C}^{M\times M}$ is the identity matrix and + denotes the
Moore–Penrose inverse.

Finally, we obtain the extracted target speech using the multichannel Wiener
filter as

$\displaystyle\hat{s}_{ij}^{(\mathrm{t})}=r_{ij}^{(\mathrm{t})}R_{i}^{(\mathrm{t})}\quantity(R_{ij}^{(\mathrm{o})})^{-1}\bm{x}_{ij}.$
(26)

##  3 Proposed Method

###  3.1 Motivation and Framework

In recent years, research on DNNs has become more popular and DNN-based source
separation methods have been proposed [19, 18, 20, 21, 17]. IDLMA [17] has
been proposed for determined audio source separation. IDLMA has higher
separation performance than ILRMA by introducing the supervised DNNs into the
source models. Since IDLMA is also a rank-1 method similarly to ILRMA, we
propose the use of IDLMA instead of ILRMA as the preprocessing method for
RCSCME. We expect that RCSCME will achieve high separation performance using
highly accurate parameters pre-estimated by IDLMA. However, since the original
IDLMA dealt with musical source separation, we should prepare new DNNs for the
speech and the noise. For the single-channel speech enhancement DNN, Denoiser
[18] has been proposed. Here, we propose a scheme to estimate the power
spectrogram not only of the target speech, but also of the noise using
Denoiser.

As shown in Table 1, Denoiser tends to have very high accuracy in terms of the
source-to-inference ratio (SIR) [22] improvement, which is a metric of
denoising, but poor results in terms of the sources-to-artifact ratio (SAR)
[22], which is a metric of artifacts. In other words, Denoiser is very good at
removing noise, but it also tends to partially remove the target speech.
Considering the high SIR of Denoiser, we can extract the noise-only time
intervals from the observed signal corresponding to the quiet time intervals
of the output of Denoiser and approximate the noise SCM from these intervals.
Thus, we propose noise self-supervised RCSCME, in which the prior distribution
of the noise SCM is introduced using the approximated noise SCM.

Table 1: SIR improvement and SAR of Denoiser | Noise
---|---
| Babble | Cafe | Station | Traffic
SIR improvement [dB] | 17.45 | 19.23 | 18.33 | 20.64
SAR [dB] | 8.89 | 9.94 | 7.54 | 10.83

Finally, an overview of the proposed method is as follows:

  1. 1.

Prepare the trained Denoiser.

  2. 2.

Estimate the demixing matrix $W_{i}$ by IDLMA using Denoiser.

  3. 3.

Calculate the quiet time intervals of the output using Denoiser, and set the
corresponding time intervals of the observed signal as the noise-only signal.

  4. 4.

Determine the prior distribution of the noise from the noise-only signal
calculated in 3) and perform RCSCME with the result of 2).

  5. 5.

Construct the multichannel Wiener filter with the result of 4).

The process flow of the proposed method is shown in Fig. 1.

Figure 1: Process flow of noise self-supervised RCSCME.

###  3.2 IDLMA for Speech Enhancement

In the separated signals $Y_{n}$ obtained by IDLMA, the $n_{t}$th signal must
be the target speech, and the other channels must be noise signals. Therefore,
for the variance of the $n_{t}$th separated signal, we use the output of the
DNN, and for the other channels, we subtract the outputs of the DNN from the
separated signals to simulate the estimation of the noise signal:

$\displaystyle\zeta_{ij,n}$
$\displaystyle=\begin{cases}\quantity(\mathrm{DNN}\quantity(Y_{n}))_{ij}&\text{if
}n=n_{t},\\\
y_{ij,n}-\quantity(\mathrm{DNN}\quantity(Y_{n}))_{ij}&\text{otherwise,}\end{cases}$
(27) $\displaystyle\sigma_{ij,n}^{2}$
$\displaystyle=\max\quantity{\quantity|\zeta_{ij,n}|^{2},\varepsilon}.$ (28)

The reason for using the DNN for each channel is that each separated signal
and each output of the DNN need to be in phase when we obtain the noise signal
by subtraction.

###  3.3 Self-Estimation of Prior Distribution of Noise

We consider the silent parts of the output of Denoiser to be the noise-only
parts in the observed signal. By thresholding the denoised signal by Denoiser,
we extract the noise-only signal
$\quantity{{\bm{x}}_{ij^{\prime}}}_{j^{\prime}\in\mathcal{J}^{\prime}}$, where

$\displaystyle\mathcal{J}^{\prime}$
$\displaystyle=\left\\{j\in\left\\{1,\cdots,J\right\\}\middle|\sqrt{\sum_{i}\absolutevalue{\quantity(\mathrm{DNN}\quantity(X_{n_{t}}))_{ij}}^{2}}<\theta\right\\}.$
(29)

Here, $\theta>0$ is a threshold parameter. The SCM of the noise-only time
intervals is calculated as

$\displaystyle\breve{R}_{i}^{(\mathrm{n})}$
$\displaystyle=\mathbb{E}\,\quantity[{\bm{x}}_{ij^{\prime}}{\bm{x}}_{ij^{\prime}}^{\mathsf{H}}]$
$\displaystyle\approx\frac{1}{\\#\mathcal{J^{\prime}}}\sum_{j^{\prime}}{\bm{x}}_{ij^{\prime}}{\bm{x}}_{ij^{\prime}}^{\mathsf{H}},$
(30)

where $\\#$ denotes the cardinality of a finite set.

One of the conjugate prior distributions for the covariance matrix in the
multivariate complex Gaussian distribution is the complex inverse matrix gamma
distribution (CIMGD) [23], which is an extension of the inverse gamma
distribution to a positive definite Hermitian matrix. We assume that
$R_{i}^{(\mathrm{n})}$ is generated by the $M$-dimensional CIMGD:

$\displaystyle
p\quantity(R_{i}^{(\mathrm{n})};\alpha^{\prime},\beta^{\prime},\breve{R}_{i}^{(\mathrm{n})})$
$\displaystyle\propto\quantity(\det
R_{i}^{(\mathrm{n})})^{-(\alpha^{\prime}+M)}$
$\displaystyle\exp(-\frac{1}{\beta^{\prime}}\tr(\breve{R}_{i}^{(\mathrm{n})}\quantity(R_{i}^{(\mathrm{n})})^{-1})),$
(31)

where $\alpha^{\prime}>M-1$, $\beta^{\prime}>0$, and
$\breve{R}_{i}^{(\mathrm{n})}$ are the shape parameter, the scale parameter,
and the scale matrix, respectively. If $\alpha^{\prime}=\tau$ and
$\beta^{\prime}=1$ in (31), this distribution coincides with the complex
inverse Wishart distribution with degree of freedom $\tau$.

###  3.4 Log-Posterior with Prior Distribution of Noise

The log-posterior likelihood given the prior distributions for the variance of
target speech and the noise SCM is

$\displaystyle\mathcal{L}(\Theta)=$
$\displaystyle\log\prod_{i,j}p(\bm{x}_{ij}|\Theta)p(r_{ij}^{(\mathrm{t})};\alpha,\beta)p\quantity(R_{i}^{(\mathrm{n})};\alpha^{\prime},\beta^{\prime},\breve{R}_{i}^{(\mathrm{n})})$
$\displaystyle=$
$\displaystyle-\sum_{i,j}\Bigg{(}\bm{x}_{ij}^{\mathsf{H}}R_{ij}^{(\mathrm{o})}\bm{x}_{ij}+\log\det
R_{ij}^{(\mathrm{o})}$ $\displaystyle+(\alpha+1)\log
r_{ij}^{(\mathrm{t})}+\frac{\beta}{r_{ij}^{(\mathrm{t})}}\Bigg{)}$
$\displaystyle-\sum_{i}\Bigg{(}(\alpha^{\prime}+M)\log\det
R_{i}^{(\mathrm{n})}$
$\displaystyle+\frac{1}{\beta^{\prime}}\tr(\breve{R}_{i}^{(\mathrm{n})}\quantity(R_{i}^{(\mathrm{n})})^{-1})\Bigg{)}+\mathrm{const.}$
(32)

We calculate the $Q$ function in the same manner as in [15].

$\displaystyle Q\quantity(\Theta;\tilde{\Theta})=$
$\displaystyle-\sum_{i,j}\Bigg{(}(\alpha+2)\log
r_{ij}^{(\mathrm{t})}+\frac{\hat{r}_{ij}^{(\mathrm{t})}+\beta}{r_{ij}^{(\mathrm{t})}}$
$\displaystyle+M\log r_{ij}^{(\mathrm{n})}+\log\det R_{i}^{(\mathrm{n})}$
$\displaystyle+\frac{\tr(\hat{R}_{ij}^{(\mathrm{n})}\quantity(R_{i}^{(\mathrm{n})})^{-1})}{r_{ij}^{(\mathrm{n})}}\Bigg{)}$
$\displaystyle-\sum_{i}\Bigg{(}\quantity(\alpha^{\prime}+M)\log\det
R_{i}^{(\mathrm{n})}$
$\displaystyle+\frac{1}{\beta^{\prime}}\tr(\breve{R}_{i}^{(\mathrm{n})}\quantity(R_{i}^{(\mathrm{n})})^{-1})\Bigg{)}+\mathrm{const.}$
(33)

It is easily seen that we only need to change the update rule of $\lambda_{i}$
in (21). Using (24) and (25), we focus on the term depending on $\lambda_{i}$
and modify (33) to

$\displaystyle Q\quantity(\Theta;\tilde{\Theta})=$
$\displaystyle\quantity(\alpha^{\prime}+M+J)\log\lambda_{i}$
$\displaystyle+\frac{1}{\lambda_{i}}\bm{u}_{i}^{\mathsf{H}}\quantity(\frac{1}{\beta^{\prime}}\breve{R}_{i}^{(\mathrm{n})}+\sum_{j}\frac{1}{{r}_{ij}^{(\mathrm{n})}}\hat{R}_{ij}^{(\mathrm{n})})\bm{u}_{i}$
$\displaystyle+\mathrm{const.}$ (34)

We obtain the following update rule, instead of (21), by solving ${\partial
Q}/{\partial\lambda_{i}}=0$:

$\displaystyle\lambda_{i}\leftarrow\frac{1}{\alpha^{\prime}+M+J}\bm{u}_{i}^{\mathsf{H}}\quantity(\frac{1}{\beta^{\prime}}\breve{R}_{i}^{(\mathrm{n})}+\sum_{j}\frac{1}{\tilde{r}_{ij}^{(\mathrm{n})}}\hat{R}_{ij}^{(\mathrm{n})})\bm{u}_{i}.$
(35)

Figure 2: Impulse response recording conditions in [16]. Figure 3: Largest SDR
improvement among the iterations under each noise. Figure 4: Behavior of SDR
improvements under each noise.

##  4 Experiment

###  4.1 Experimental Conditions

We conducted a simulated experiment to verify the effectiveness of the
proposed method. The diffuse noise was simulated by convolution of the impulse
response from 19 loudspeakers to four microphones. We used the four noise
signals: babble, cafe, station, and traffic noise. For the babble noise, we
used 19 speech signals from the JNAS speech corpus [24]. For the cafe,
station, and traffic noise, we used 19 signals, which were obtained by
dividing each noise signal from DEMAND [25]. The target speech signal was
simulated by convolution of the impulse response from a closer location. For
the target speech signal, we used one speech signal from the JNAS speech
corpus [24]. We used the conditions in Fig. 2 to record the impulse response
[16]. The sampling rate was 16 kHz. The input signal-to-noise ratio was set to
0 dB. The STFT was performed by using a 64-ms-long Hamming window and 32-ms-
long shift.

We compared ILRMA, IDLMA, RCSCME using ILRMA (ILRMA+RCSCME), RCSCME using
IDLMA with $\sigma_{ij,n}$ updated by (27) and (28) (proposed method 1), and
noise self-supervised RCSCME using IDLMA with $\lambda_{i}$ updated by (35)
(proposed method 2). We used source-to-distortion ratio (SDR) [22] improvement
as a measure. We used
$(\alpha,\beta,\alpha^{\prime},\beta^{\prime})=(1.3,10^{-16},8\times
10^{2},10^{4})$ for each method, which were selected experimentally. The
flooring parameter in (6) and (28) was set to
$\varepsilon=0.1\times\sum_{i,j}\sigma_{ij,n}^{2}/(IJ)$ and the threshold
parameter in (29) was set to $\theta=10^{-3}$. The demixing matrix $W_{i}$ in
ILRMA and IDLMA was initialized by the identity matrix. The number of NMF
bases of ILRMA was set to $10$ and the NMF variables were initialized by
random values taken from the uniform distribution on $[0,1]$. We used the unit
eigenvector corresponding to the zero eigenvalue of
${{R^{\prime}}_{i}^{(\mathrm{n})}}$ as $\bm{v}_{i}$ in RCSCME. For Denoiser,
we used the trained model provided by [18] with $H=48$ hidden channels. Since
Denoiser is the waveform-domain model, we used Denoiser as

$\displaystyle\mathrm{DNN}(Y_{n})=\mathrm{STFT}\quantity[\mathrm{Denoiser}\quantity(\mathrm{STFT}^{-1}\quantity[Y_{n}])].$
(36)

We set the number of iterations in ILRMA to 50, that in IDLMA to 90, and that
in RCSCME for ILRMA+RCSCME, proposed method 1, and proposed method 2 to 10. In
IDLMA, $\sigma_{ij,n}$ was updated by DNN for every 30 times $W_{i}$ was
updated. These numbers of iterations were experimentally determined.

###  4.2 Results

We took the average of the SDR improvement results for 10 different sets of
randomly initialized values of the NMF variables in ILRMA and ILRMA+RCSCME.
The SDR improvements for each noise are shown in Fig. 3. For each RCSCME
method, we show the largest SDR improvement among the iterations. As shown in
Fig 3, Denoiser outperformed IDLMA for babble, cafe, and traffic noise. The
reason is that Denoiser includes a nonlinear transformation that can provide
stronger noise reduction than linear systems, whereas IDLMA estimates the
output signals via a linear time-invariant spatial filter, i.e., demixing
matrix $W_{i}$. However, proposed method 1 shows better performance than
Denoiser and IDLMA for all the noises, implying the advantage of the time-
variant properties in RCSCME. It is also shown that proposed method 2 achieves
the best separation performance among all methods; this would be due to
introduction of the self-estimated prior distribution of the noise SCM. The
SDR improvement for each iteration is shown in Fig. 4. Fig. 4 shows that
ILRMA+RCSCME requires more than five iterations to converge, while proposed
method 1 and proposed method 2 take the maximum value in less than five
iterations.

##  5 Conclusions

We proposed RCSCME using IDLMA, which is applied to speech enhancement by
preparing a single-channel speech enhancement DNN for estimating the power
spectrogram of the target speech and the noise. We used Denoiser as the
single-channel speech enhancement DNN model for estimating the power
spectrogram of the target speech. To estimate the power spectrogram of diffuse
noise, we subtracted the speech extracted by Denoiser from the separated
signal. We also proposed a method to estimate the noise-only time intervals
using the output of Denoiser and design the prior distribution of the noise
SCM for RCSCME. We confirmed that noise self-supervised RCSCME using IDLMA
outperformed the other methods and took the maximum value in fewer iterations
than the conventional methods.

## Acknowledgment

This work was supported by the Japan–New Zealand Research Cooperative Program
of JSPS and RSNZ (Grant Number JPJSBP120201002), JSPS KAKENHI Grant Numbers
19K20306, 19H01116, and 19H04131, and JST Moonshot R&D Grant Number
JPMJPS2011.

## References

  * [1] H. Sawada, N. Ono, H. Kameoka, D. Kitamura, and H. Saruwatari, “A review of blind source separation methods: two converging routes to ILRMA originating from ICA and NMF,” _APSIPA Trans. Signal Inf. Process._ , vol. 8, no. e12, pp. 1–14, 2019.
  * [2] A. Hiroe, “Solution of permutation problem in frequency domain ICA, using multivariate probability density functions,” in _Proc. ICA_ , 2006, pp. 601–608.
  * [3] T. Kim, H. T. Attias, S.-Y. Lee, and T.-W. Lee, “Blind source separation exploiting higher-order frequency dependencies,” _IEEE Trans. ASLP_ , vol. 15, no. 1, pp. 70–79, 2006.
  * [4] N. Ono, “Stable and fast update rules for independent vector analysis based on auxiliary function technique,” in _Proc. WASPAA_ , 2011, pp. 189–192.
  * [5] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and H. Saruwatari, “Determined blind source separation unifying independent vector analysis and nonnegative matrix factorization,” _IEEE/ACM Trans. ASLP_ , vol. 24, no. 9, pp. 1626–1641, 2016.
  * [6] D. D. Lee and H. S. Seung, “Learning the parts of objects by non-negative matrix factorization,” _Nature_ , vol. 401, no. 6755, pp. 788–791, 1999\.
  * [7] S. Araki, R. Mukai, S. Makino, T. Nishikawa, and H. Saruwatari, “The fundamental limitation of frequency domain blind source separation for convolutive mixtures of speech.” _IEEE Trans. SAP_ , vol. 11, no. 2, pp. 109–116, 3 2003.
  * [8] S. Araki, S. Makino, Y. Hinamoto, R. Mukai, T. Nishikawa, and H. Saruwatari, “Equivalence between frequency-domain blind source separation and frequency-domain adaptive beamforming for convolutive mixtures,” _EURASIP J. Applied Signal Process._ , vol. 2003, no. 11, pp. 1157–1166, 2003\.
  * [9] Y. Takahashi, T. Takatani, K. Osako, H. Saruwatari, and K. Shikano, “Blind spatial subtraction array for speech enhancement in noisy environment.” _IEEE Trans. ASLP_ , vol. 17, no. 4, pp. 650–664, 5 2009.
  * [10] A. Ozerov and C. Févotte, “Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation,” _IEEE Trans. ASLP_ , vol. 18, no. 3, pp. 550–563, 2009.
  * [11] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Multichannel extensions of non-negative matrix factorization with complex-valued data,” _IEEE Trans. ASLP_ , vol. 21, no. 5, pp. 971–982, 2013.
  * [12] N. Ito and T. Nakatani, “FastMNMF: Joint diagonalization based accelerated algorithms for multichannel nonnegative matrix factorization,” in _Proc. ICASSP_ , 2019, pp. 371–375.
  * [13] K. Sekiguchi, A. A. Nugraha, Y. Bando, and K. Yoshii, “Fast multichannel source separation based on jointly diagonalizable spatial covariance matrices,” in _Proc. EUSIPCO_ , 2019, 5 pages.
  * [14] N. Q. Duong, E. Vincent, and R. Gribonval, “Under-determined reverberant audio source separation using a full-rank spatial covariance model,” _IEEE Trans. ASLP_ , vol. 18, no. 7, pp. 1830–1840, 2010.
  * [15] Y. Kubo, N. Takamune, D. Kitamura, and H. Saruwatari, “Efficient full-rank spatial covariance estimation using independent low-rank matrix analysis for blind source separation,” in _Proc. EUSIPCO_ , 2019, pp. 1–5.
  * [16] Y. Kubo, N. Takamune, D. Kitamura, and H. Saruwatari, “Blind speech extraction based on rank-constrained spatial covariance matrix estimation with multivariate generalized Gaussian distribution,” _IEEE/ACM Trans. ASLP_ , vol. 28, pp. 1948–1963, 2020.
  * [17] N. Makishima, S. Mogami, N. Takamune, D. Kitamura, H. Sumino, S. Takamichi, H. Saruwatari, and N. Ono, “Independent deeply learned matrix analysis for determined audio source separation,” _IEEE/ACM Trans. ASLP_ , vol. 27, no. 10, pp. 1601–1615, 2019.
  * [18] A. Defossez, G. Synnaeve, and Y. Adi, “Real time speech enhancement in the waveform domain,” in _Proc. Interspeech_ , 2020.
  * [19] E. M. Grais, M. U. Sen, and H. Erdogan, “Deep neural networks for single channel source separation,” in _Proc. ICASSP_ , 2014, pp. 3734–3738.
  * [20] K. Qian, Y. Zhang, S. Chang, X. Yang, D. Florencio, and M. Hasegawa-Johnson, “Deep learning based speech beamforming,” in _Proc. ICASSP_ , 2018, pp. 5389–5393.
  * [21] A. A. Nugraha, A. Liutkus, and E. Vincent, “Multichannel audio source separation with deep neural networks,” _IEEE/ACM Trans. ASLP_ , vol. 24, no. 9, pp. 1652–1664, 2016.
  * [22] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” _IEEE Trans. ASLP_ , vol. 14, no. 4, pp. 1462–1469, 2006.
  * [23] A. Iranmanesh, M. Arashi, D. Nagar, and S. Tabatabaey, “On inverted matrix variate gamma distribution,” _Commun. Stat. - Theory Methods_ , vol. 42, no. 1, pp. 28–41, 2013.
  * [24] K. Itou, M. Yamamoto, K. Takeda, T. Takezawa, T. Matsuoka, T. Kobayashi, K. Shikano, and S. Itahashi, “JNAS: Japanese speech corpus for large vocabulary continuous speech recognition research,” _JASJ (E)_ , vol. 20, no. 3, pp. 199–206, 1999.
  * [25] J. Thiemann, N. Ito, and E. Vincent, “DEMAND: a collection of multi-channel recordings of acoustic noise in diverse environments,” Zenodo, Jun. 2013. [Online]. Available: https://doi.org/10.5281/zenodo.1227121