ICASSP 2018, Calgary, Canada

Tutorial T1(April 16th, 2018)

Blind Audio Source Separation
on Tensor Representation
Hiroshi Sawada (NTT Corporation)
Nobutaka Ono (Tokyo Metropolitan University)
Hirokazu Kameoka (NTT Corporation)
Daichi Kitamura (The University of Tokyo)
1

Tutorial structure
1. Introduction
1. Separation of audio/speech signals
2. Live demonstration
2. ICA and IVA
1. ICA: Independent Component Analysis
2. IVA: Independent Vector Analysis
3. NMF
1. NMF: Nonnegative Matrix Factorization
2. MNMF: Multichannel NMF
4. ILRMA
1. ILRMA: Independent Low-Rank Matrix Analysis

2Separation of audio/speech signal
cocktail party effectspeech recognition
in noisy environment

music analysis

3Live demonstration• Separate 2 speeches with 2 microphones• iPhone app• Script
The ICASSP meeting is theworld’s largest and mostcomprehensive technicalconference focused on signal
processing and its applications.

We are demonstrating theBlind Source Separation forconvolutive mixtures ofspeech in a real room withreal talkers.

4

BSS: Blind Source Separation• Separate the mixtures at microphones• into the original sources• assuming M=N=2 for simplicity• source activity and mixing system H is unknown (blind)
+

+

+

+

5

Instantaneous BSS• Mixing system H is described by scalars• Sources are multiplied by scalars and then mixed+
Simple and basic
BSS model+

6

Convolutive BSS• Delay and reverberations in a real room situation• mixing system H is described by impulse responses• Sources are convolutively mixed+
Convolutive BSS is a muchharder problem thaninstantaneous BSS

+

7

8

The whole system
Sources

Mixtures

Separations

+

+

+

+

Mixing system

Separation system

STFT: short-time Fourier transformconvolution

Convolution filter(time domain)
STFT

Frequency response(frequency domain)

multiplication

STFT: short-time Fourier transform• From a time-domain real-valued signal• To a time-frequency-domain complex-valued signal
STFT
Spectrogram: only amplitudes are displayed

910

Three key axes
Frequency• Frequency-domain processing is effective• source characteristic• convolution → multiplication
Time• Source activity, Onset and offset
ChannelChannel• Source, Mixture, Separation
MixtureTensor representation Frequency
Time

11

Tensor and sliced matrices
2.2. IVA
4.1. ILRMA

Channel
2.1. ICA
Mixturevector ← component

ChannelMixtureTime

slice into matrices
Frequency
Time
Frequency

3.2. MNMF
multichannel

3.1. NMF

TimeNotations
s: sources
x: mixtures/observations
y: separations
W: separation system
H: mixing system

i: frequency bin index
j: time frame index
m: microphone index
n: source/separation index
I: number of frequency bins
J: number of time frames
M: number of microphones
N: number of sources

T: basis spectrum
V: time-varying magnitude (NMF), whitening matrix (FastICA),
weighted covariance matrix (AuxIVA)

12

Tutorial structure
1. Introduction
1. Separation of audio/speech signals
2. Live demonstration
2. ICA and IVA
1. ICA: Independent Component Analysis
2. IVA: Independent Vector Analysis
3. NMF
1. NMF: Nonnegative Matrix Factorization
2. MNMF: Multichannel NMF
4. ILRMA
1. ILRMA: Independent Low-Rank Matrix Analysis

13

14

Tensor and sliced matrices
Channel
2.1. ICA
MixtureChannelMixtureTime

slice into matrices
Frequency
Time

Convolutive BSS

becomes easier

Instantaneous BSS

permutation indeterminacy needs to be aligned
to properly reconstruct the whole separations

ICA: Independent Component Analysis• Extract the original sources from the mixtures• Mixing matrix H cannot be obtained
+

+

+

+

Sources are
assumed to be
independent

Makes
independent toeach other15

ICA: Independent Component Analysis• In addition to independence, need to assume• source distributions are different from Gaussian
+

+

+

+

Distributions aredifferent fromGaussian

Makes thedistribution far
from Gaussian

16

17

Speech mixtures
amplitudessamplitudeamplitude

18sStatistics of the mixtures

amplitude

19

Statistics of the mixtures

amplitude→ converge to the red line (Gaussian)Central Limit Theorem20

A complex-valued source model• Super-Gaussian distribution• Sharper peak at the origin than Gaussian

Real
part

Imaginarypart21Independent component analysis+

Unknown
mixingprocess+• Linear operation• Output independence• Non-gaussianity

22Likelihood of separation matrix W• Likelihood of W for the whole observations• Probability density function, linear transformation• Output independence
Log-likelihood function

23

Objective function to be minimized

24

Maximum likelihood estimation = Minimize the negative log-likelihood

Contrast function

1st order derivative

2nd order derivativeOptimization methods• Gradient descent:step size

• Expensive matrix inversion• Slow convergence• Three practical ways
1. Natural gradient2. Pre-whitening + FastICA
3. Auxiliary function-based optimization

25

26

1. Natural gradient[Amari et al., 1996]• No matrix inversion• Efficient computation• Equivariance property [Cardoso and Souloumiac, 1996]• Free from the characteristics of mixing matrix(e.g. close to singular)

Convergence example
Estimated mixing matrix
Starting from

Step size:27

2. Pre-whitening ＋ FastICA[Hyvarinen et al., 2001]

Separation matrix of the form:
Mixtures

Pre-whitenings. t.
via eigenvalue
decomposition of

28

Whitened

Separated

FastICA
Minimize the objective function
w.r.t. a unitary matrix U

FastICA• Objective function w.r.t a unitary matrix U• Minimizefor• with unitary constraint

29

FastICA algorithm

30[Hyvarinen et al., 2001]• For(sequentially)• Iterate the followings until convergence

Separated signal calculation

Optimization of G by Newton’s method

Gram-schmidt orthogonalization
Unit-norm normalization

FastICA convergence example• Red（□）• Starting from
• Likelihood maximization: points close to the origin• Unit-norm normalization: points on the unit sphere• Good solution only by 5 iterations• Green（△）• Starting from• One-step solution by orthogonalization

31

3. Auxiliary function• Objective function to be minimized
• If

and
is monotonically decreasing• Auxiliary function
with auxiliary variable

Equal when

32

[Ono and Miyabe, 2010][Ono 2011]

Auxiliary function-based optimization• Iterate the followings until convergence
Separated signal calculation
Weighted covariance matrices for each separation

Solve the HEAD problem for33[Ono and Miyabe, 2010][Ono 2011]

Auxiliary ICA convergence example
Estimated mixing matrix
Starting fromable to take a big stepat an early stage

34

Comparisons of three methods

35• measured by signal-to-interference ratio (SIR)• at frequency 3586 Hz, 201 samples, pre-whitening applied

convergence behavior

computational cost considered Natural gradient: sensitive to step-size FastICA: fast convergence, slightly limited SIR (unitary constraint) Auxiliary ICA: fast convergencePermutation and Scaling problem• Ambiguities of ICA solutions
If

is a solution, thenis also a solution

for any diagonal

Independence ofand permutation

matrixdoes not change

36Permutation and Scaling problem

Frequency bin

Raw ICA outputs

Time frame

PermutationalignedPermutation/Scalingaligned37

Solving Permutation and Scaling Problems

38• Permutation• Post-processing• see e.g. [Sawada et al., 2004], [Sawada et al., 2011]• Tensor methods (IVA, ILRMA)• will be explained in later sections• Scaling• Refer to a microphone observation[Cardoso 1998][Murata et al., 2001]• So-called “projection back”
Nakashima 2001]• via mixing system estimation [Matsuoka and[Takatani et al., 2004]

39

Scaling alignment example
Untouched

Alignedreferencemicrophone

estimated mixing scalarMixing system estimation
Estimated mixing situation

ICA result

How to calculate matrix
If

has an inverse

Otherwise ()• Least-mean-square estimatorthat minimizes• Moore-Penrose pseudo inverse

40

Tutorial structure
1. Introduction
1. Separation of audio/speech signals
2. Live demonstration
2. ICA and IVA
1. ICA: Independent Component Analysis
2. IVA: Independent Vector Analysis
3. NMF
1. NMF: Nonnegative Matrix Factorization
2. MNMF: Multichannel NMF
4. ILRMA
1. ILRMA: Independent Low-Rank Matrix Analysis

41

42

Tensor and sliced matrices
2.2. IVA

Channel
2.1. ICA
Mixturevector ← component

ChannelMixtureFrequency
Time

Time

Two papers at ICA 2006

Concept of “multivariate source model” was presented
in two papers at the same conference independently

43Independent vector analysis• ICA: Sources generate stochastic scalar variables
Source 1
Source 2

MixtureChannelFrequency

Time

Mixing process:Separation process:• IVA: Sources generate stochastic vector variables
Source 1
Source 2

MixtureSame mixing / separation process as ICA, but multivariate source model

44

Likelihood of separation matrices W

45• Likelihood of all Ws for the whole observations• Probability density function, linear transformation• Output independence
ChannelMixtureLog-likelihood function

Frequency
TimeObjective function of IVA• A set of demixing matrices to be estimated• Objective function of IVA

Contrastfunction
Multivariate p.d.f.

FrequencyTimey is a function of w

46

What p.d.f. is appropriate for a spectrum vector?

47

Choice of multivariate p.d.f.
• Necessary properties• Non-Gaussian (like ICA)
• Representing higher-order dependency
between vector components
FrequencyTime• Well-used multivariate p.d.f.
• Spherical super-Gaussian [Hiroe 2006],[Kim 2006]• Time-varying Gaussian [Ono+ 2012]
Variance is time-varying.
Totally it is super-Gaussian.
(Show later)

48

49

Why spherical super-Gaussian?
• Let

where

• Whenis Gaussian, sub-Gaussian, and superGaussian, what dependency between
andisrepresented?
0.6

0.2

0.8

0.15

0.6

)
2

,y

0.05

p(y

p(y

0.2

0

0

0

2
-2

y
1

-2
0

2

0

2

-2
1

2

Gaussian

0

2

0

2

y
y

0.2

0
-2

-2

0.4
1

1

p(y

0.1
1

,y

,y

2

2

)

)

0.4

-2

y

y
2

Sub-Gaussian

0

2
1

y
2

Super-Gaussian

50

Spherical Gaussian
Joint p.d.f.

Conditional p.d.f.
0.7

0.6
0.6

0.5

0.4

|y

Completely
overlapped

2

0.2

p(y

p(y

1

1

)

,y

2

)

0.4

0.3

0.2

0
0.1

-2
0

0

2
0

2

-3

-2

-1

0

1

2

y
2

-2

y
1

y
2

In spherical Gaussian case, the value of y1 does not change
the p.d.f. of y2. It means y1 and y2 do not have any dependencies.

3

51

Spherical Sub-Gaussian
Joint p.d.f.

Conditional p.d.f.
0.6

0.2
0.5

0.15

2

)

0.4

)
|y
2

0.3

0.05

p(y

p(y

1

1

,y

0.1

0.2

0
0.1

-2
0

0

2
0

2

-3

-2

-1

0

1

y
2

-2

y
1

y
2

In spherical sub-Gaussian case, when |y1| is larger,|y2| tends to be smaller oppositely.

2

3

52

Spherical Super-Gaussian
Joint p.d.f.

Conditional p.d.f.
1

0.8

0.8

0.4

0.6

)
1
2

0.2

p(y

p(y

|y

1

,y

2

)

0.6

0.4

0
0.2

-2
0

0

2
0

2
1

-2

-1

0

1

y
2

-2

y

-3

y
2

In spherical super-Gaussian case, when |y1| is larger,|y2| tends to be also larger. Therefore, it representsco-occurrence among components.

2

3

53

Time-Varying GaussianFrequency• P.d.f. of time-varying Gaussian• Co-occurrence among frequency
components are explicitly
represented by.
• Shared among all frequency i.
• Can be changed at each time frame j.

TimeSolutions for IVA

54• The objective function of IVA is also nonlinear.

• Similarly as ICA, three typical methods• Natural gradient [Kim+ 2006,Hiroe 2006]• Pre-whitening + Fixed-point iteration (FastIVA) [Lee+ 2007]• Dr. Taesu Kim (an inventor of IVA)’s code is available
https://github.com/teradepth/iva• Auxiliary function-based optimization (AuxIVA) [Ono2011, Ono2012b]

Auxiliary Function Approach• Optimization Problem• Auxiliary Function• Alternative Update Rules
Auxiliary variable update (like E-step)
Parameter update (like M-step)• Advantages• Stable: Convergence is guaranteed• Simple: No tuning parameters such as step size
But how to find useful auxiliary function is problem-dependent

More details will be explained later

55

Theorem for quadratic auxiliary function
• If

and

is monotonically decreasing inSpherical,

Super-Gaussian
Nonlinear
contrast function

Quadratic
auxiliary function

holds. The equality sign is valid iff
Examples56

Multivariate
p.d.f.

Contrastfunction.
Weightfunction

0

r

0

r

0

r

Auxiliary Function for IVA• Objective Function of IVA

If G is spherical andderived from super-Gaussian

Nonlinearfunction of w

• Auxiliary Function for IVA
Quadratic
function of w

Source activity(auxiliary variable)
Weightedcovariance matrix:

57

Minimizing auxiliary function• The demixing matrix should be updated such that auxiliary
function is minimized.

• Fromfor all,

HEAD (Hybrid Exact-Approximate Joint Diagonalization)
problem [Yeredor 2009] is derived.

58

HEAD Problem (1/2)

59

For simplicity, a frequency index i is dropped in this slide.

Given N positive definite matrices
find an N×N matrix,
such that[Yeredor 2009]

ex. N=3 case60

HEAD Problem (2/2)
For simplicity, a frequency index i is dropped in this slide.

Given N positive definite matrices
find an N×N matrix,
such that[Yeredor 2009]• Remarks• Number of equations = number of variables• When, it is equivalent to generalized eigenvalue problem,
which can be solved in a closed form. (Show later)• When, a closed-form solution has never been found.

Derivation of sequential update rule61

For simplicity, a frequency index i is dropped in this slide.

Let’s consider to solve

in with fixing

Some vector linearly
independent of w2 and w3.

If we choose
obtainedat the previous iteration as
can be rewritten by.

,

We can update
sequentially in the same way.

Scale can be fixed by

62

Algorithm of AuxIVA
Iterate until convergence
For[Ono 2011](every source)
Update of separation
Update of source activity(shared in all frequency)For(every frequency)
Update ofweightedcovariance matrix
Update ofdemixing matrix
Unit vector with the nth element unity

Separation Example (1/6)

63Separation Example (2/6)

64Separation Example (3/6)

65Separation Example (4/6)

66Separation Example (5/6)

67Separation Example (6/6)

68Comparison with Natural GradientThree sources case
One iteration(Matlab)AuxIVA : 0.34s
GradIVA: 0.16s

69

HEAD Problem in Two-Sources Case

HEAD problem is deformed to generalized eigenvalue problem,
which can be solved in a closed form [Yoshioka 2008,Ono 2010]

70

71

Algorithm of Stereo AuxIVA
Iterate until convergence[Ono 2012b]

Update of separation
Update of source activity(shared in all frequency)

Update of weighted
covariance matrix
Find

by solving

Update of demixing matrix

Comparison with Natural Gradient (2)Two sources case

SIR improvement [dB]15

10Natural Grad. (mu=0.1)
Natural Grad. (mu=0.2)
Natural Grad. (mu=0.3)
AuxIVA1
Stereo AuxIVA
AuxIVA2

5

0
0

50

100
Iteration steps

150

200

72

73

Explicit Representation of Eigenvectors of 2×2 Matrix
For simplicity, a frequency index i is dropped in this slide.

•• Two eigenvectors of 2×2 matrix can be explicitly given byalways holds• See [Ono 2012b] about fast implementation using vector
operation.

Implementation on iPhone• Stereo AuxIVA was implemented on iPhone [Ono 2012b]• Calculation time is almost
linear to input signal length(RTF≒1/5 @ 16kHz, 10itr.
on iPhone4)• Demo on youtube(https://www.youtube.com/watch?v=lLMbflDMMeE)• Another Demo in Demo Session II at this ICASSP• DEMO-2.2: Low-Latency Real-Time Blind Source Separationfor Hearing Aids
• Time: Wednesday, April 18, 13:30 - 15:30• Location: Exhibit Hall Foyer

74

Tutorial structure
1. Introduction
1. Separation of audio/speech signals
2. Live demonstration
2. ICA and IVA
1. ICA: Independent Component Analysis
2. IVA: Independent Vector Analysis
3. NMF
1. NMF: Nonnegative Matrix Factorization
2. MNMF: Multichannel NMF
4. ILRMA
1. ILRMA: Independent Low-Rank Matrix Analysis

75

76

Tensor and sliced matrices

ChannelMixtureFrequency
Time
Frequency

3.1. NMF

TimeFrequency

What is NMF? From “generative model” perspective

time

77

What is NMF? From “generative model” perspective
Time-varying amplitude of each basis spectrum

Frequency

Basis spectrum

time

78

What is NMF? From “generative model” perspective

Frequency

Matrix factorization of 𝐗𝐗 = Restoring information about 𝐓𝐓 and 𝐕𝐕 from 𝐗𝐗

20

40

60

80

time 120

100

79

80

NMF as spectrogram model fitting• Model a mixture spectrum as the sum of basis spectra
scaled by time-varying magnitudes
Time-varying magnitude
Basis spectrum

Frequency
Time

Matrix notation:

Time-varying
magnitudes

𝑘𝑘Model-fitting
problem:Defining
divergencemeasure leads
to optimization
problemFrequency 𝑖𝑖

Basis spectra

𝑘𝑘time 𝑗𝑗

Popular divergence measures• Measure of difference between 𝑥𝑥 and 𝑦𝑦

Euclidean distance
Kullback-Leibler (KL)divergence

Itakura-Saito (IS)divegence81

Geometrical understanding of NMF• Because of the non-negativity of 𝐓𝐓, all basis vectors lie in the first quadrant.
• Because of the non-negativity of 𝐕𝐕, 𝐓𝐓𝐯𝐯𝑗𝑗 can only cover the area (a convexcone) enclosed by the extended lines of all the basis vectors.
• NMF attempts to find a convex cone that is closest to all the observed vectors.

82

Sparsity-inducing effect of NMF• NMF naturally produces sparse representations• Let 𝐯𝐯� be the solution to an unconstrained optimization problem:• 𝐓𝐓�
𝐯𝐯 corresponds to the closest point from 𝐱𝐱 in the subspacespanned by 𝐭𝐭1 , … , 𝐭𝐭 𝐾𝐾 .
• Except for the case where 𝐯𝐯� is non-negative, the solution to the� will be the closest point to𝐓𝐓𝐯𝐯�
constrained optimization problem, 𝐓𝐓𝐯𝐯,
in the area enclosed by the extended lines of all the basis vectors.
• This means at least one of the elements of 𝐯𝐯� becomes 0.
(unconstrained)(constrained)

Subspace spannedby 𝐭𝐭1 and 𝐭𝐭 283

Itakura-Saito divergence NMF [Févotte+2009]84• Model mixture signal as the sum of Gaussian-distributed
random signals with rank-1 power spectrograms
Complex spectrogram of component 𝑘𝑘

If

and

Power spectrogram of component 𝑘𝑘

Rank-1 spectrogram

are independent
Likelihood function of 𝑻𝑻 and 𝑽𝑽• Maximum likelihood of 𝐓𝐓 and 𝐕𝐕 amounts to NMF using ItakuraSaito divergence
Itakura-Saito divergenceLocal Gaussian model (LGM)• Generative model assuming each element of a complex spectrogram
to independently follow a zero-mean complex Gaussian distribution
with a different variance (power)
variance

Likely to generatecomplex numbers near 0

Frequency 𝑖𝑖

grayscaleTime 𝑗𝑗

Likely to generate complex
numbers with larger magnitudes

85

Itakura-Saito divergence NMF (IS-NMF)• Optimization problem:
Minimizesubject to• How can we solve this?
• EM algorithm [Févotte+2009]• Auxiliary function approach–Majorization-Maximization (MM) algorithm[Kameoka+2006], [Nakano+2010], [Févotte+2011]–Majorization-Equalization (ME) algorithm [Févotte+2011]86Expectation-Maximization (EM) algorithm• When a certain likelihood function can be written as
where 𝐜𝐜 is a set of latent variables,
a stationary point of
can be found by iteratively
performing the following steps:• Expectation-step• Maximization-step

Q function

87

IS-NMF optimization with EM algorithm[Févotte+2009]• Likelihood function for IS-NMF:
where• Q function(independent of 𝜃𝜃)• E-stepwhere• M-step88

Auxiliary function approach89• Techniques to find a stationary point of an objective function
using an auxiliary functionthat satisfies• Majorization-minimization

auxiliaryobjective[Hunter & Lange 2004]

[1]
[2]

[1]
[2]

• Majorization-equalization[Févotte & Idier 2010]

[1]
[2]

Graphical illustration of MM and ME algorithms

90• Comparison between MM and ME updates
auxiliary function objective function

MMMEGraphical illustration of MM and ME algorithms

91• MM algorithm
auxiliary function
objectivefunctionGraphical illustration of MM and ME algorithms

92• ME algorithm
auxiliary function
objectivefunction

Relation to EM algorithm• Objective: maximize 𝑝𝑝(𝐱𝐱|𝜃𝜃) w.r.t. 𝜃𝜃

𝐜𝐜: latent variable

Jensen’s inequality
Auxiliary function
E step：Maximize auxiliary function w.r.t. 𝑞𝑞

M step：Maximize auxiliary function w.r.t. 𝜃𝜃93Another look at MM algorithm• Coordinate descent ofparameter 𝜃𝜃Contour line ofInitialparameter
auxiliary variable 𝛼𝛼94

Another look at MM algorithm• Coordinate descent ofparameter 𝜃𝜃Contour line ofInitialparameter
auxiliary variable 𝛼𝛼95

Another look at MM algorithm• Coordinate descent ofparameter 𝜃𝜃Contour line ofInitialparameter
auxiliary variable 𝛼𝛼96Another look at MM algorithm• Coordinate descent ofparameter 𝜃𝜃Contour line ofInitialparameter
auxiliary variable 𝛼𝛼97Motivations for auxiliary function approach• Auxiliary function can be useful and effective
when one wants to• handle a non-convex objective function with multiple
local optima and stationary points,• handle an objective function that has
discontinuous/non-differentiable points,• handle equality/inequality constraints, and• avoid matrix inversions.

98Useful inequalities for auxiliary function design (1/7)• Jensen's inequality for non-negative arguments：convex function

Equality holds whenE.g.) when:right-handside

left-handside99Useful inequalities for auxiliary function design (2/7)• Jensen's inequality for real number arguments：convex function

Equality holds when

• Used for complex NMF[Kameoka+2009]

right-handside

and RBM optimization[Kameoka+2014][Takamune+2014]

left-hand
side

100

Useful inequalities for auxiliary function design (3/7)

101• 1st order Taylor expansion of convex/concave functions: concave function
Equality holds whenE.g.) when:Useful inequalities for auxiliary function design (4/7)• 1st order Taylor expansion of logarithmic function
Scalar case:
Equality holds when

Extension to matrix case:
Equality holds whenEquals to the sum of the logarithmsof the eigenvalues of 𝐗𝐗

Used in multichannel NMF frameworks [Sawada+2012, Higuchi+2014]
and Positive Semidefinite Tensor Factorization [Yoshii+2013]102

Useful inequalities for auxiliary function design (5/7)• Quadratic function tangent to power functions
When 0 < 𝑝𝑝 ≤ 2:
Equality holds when

• Used for sparse regularization
for complex NMF [Kameoka+2009]103

Useful inequalities for auxiliary function design (6/7)

104• 1st order Taylor expansion of L2 norm
When

Equality holds when

For complex number:
When• Used for sound source localization[Ono+2009][Ono+2010]

and Time-domain Spectrogram
Factorization (TSF) [Kameoka2015]

Equality holds whenUseful inequalities for auxiliary function design (7/7)• Logistic function

Equality holds when

Gaussian distribution function105

IS-NMF optimization with MM algorithm

106[Kameoka+2006][Nakano+2010][Févotte+2011]• Objective function to be minimized:



Reciprocal function is convex in positive domain

Jensen's inequality


Logarithmic function is concave
Right-handsideIS-NMF optimization with MM algorithm

107[Kameoka+2006][Nakano+2010][Févotte+2011]• Majorizer• Update rules for 𝝀𝝀 and 𝐔𝐔• Update rules for 𝐓𝐓 and 𝐕𝐕* Update rules of ME algorithm can be obtained by solving 𝐓𝐓 and 𝐕𝐕
that satisfies

Tutorial structure
1. Introduction
1. Separation of audio/speech signals
2. Live demonstration
2. ICA and IVA
1. ICA: Independent Component Analysis
2. IVA: Independent Vector Analysis
3. NMF
1. NMF: Nonnegative Matrix Factorization
2. MNMF: Multichannel NMF
4. ILRMA
1. ILRMA: Independent Low-Rank Matrix Analysis

108

109

Tensor and sliced matrices

ChannelMixtureFrequency
Time
Frequency

3.2. MNMF
multichannel

3.1. NMF

TimeExtension to multichannel input

110• Frequency-domain instantaneous mixture:

Mixing matrix

Source signals

channel𝑚𝑚

Observed signals

frequency
𝑖𝑖

time
𝑗𝑗

source𝑛𝑛time𝑗𝑗

source𝑛𝑛• Assume local Gaussian model (LGM) with source power
spectrograms expressed using NMF (low-rank) modelLocal Gaussian model• Source model assuming each element of a complex spectrogram toindependently follow a zero-mean complex Gaussian distribution witha different variance (power)
variance

111

112

Local Gaussian model• Source model assuming each element of a complex spectrogram toindependently follow a zero-mean complex Gaussian distribution witha different variance (power)
variance

Likely to generatecomplex numbers near 0

Frequency 𝑖𝑖

grayscaleTime 𝑗𝑗

Likely to generate complex
numbers with larger magnitudes

Local Gaussian model• Source model assuming each element of a complex spectrogram toindependently follow a zero-mean complex Gaussian distribution witha different variance (power)
variance2• Allows us to incorporate power spectrogram model in 𝜎𝜎𝑖𝑖𝑖𝑖,𝑛𝑛2• 𝜎𝜎𝑖𝑖𝑖𝑖,𝑛𝑛= ∑𝑘𝑘 𝑡𝑡𝑖𝑖𝑘𝑘,𝑛𝑛 𝑣𝑣𝑘𝑘𝑗𝑗,𝑛𝑛 corresponds to NMF model

113

114

Local Gaussian model• Source model assuming each element of a complex spectrogram toindependently follow a zero-mean complex Gaussian distribution witha different variance (power)
variance

Likelihood

frequency𝑖𝑖

time 𝑗𝑗

source𝑛𝑛time 𝑗𝑗

source 𝑛𝑛channel𝑚𝑚2• Allows us to incorporate power spectrogram model in 𝜎𝜎𝑖𝑖𝑖𝑖,𝑛𝑛2• 𝜎𝜎𝑖𝑖𝑖𝑖,𝑛𝑛= ∑𝑘𝑘 𝑡𝑡𝑖𝑖𝑘𝑘,𝑛𝑛 𝑣𝑣𝑘𝑘𝑗𝑗,𝑛𝑛 corresponds to NMF model

Super-Gaussianity of local Gaussian model

115• Theorem
The time average of zero-mean Gaussian distributions
with time-varying variances is super-Gaussian.

• Proof omitted
0
Kurtosis becomes 0 if and only
if all the variances are equal.

Using local Gaussian models implies assuming
source signals to follow super-Gaussian distributions

[Ozerov+2010, Sawada+2012,...]

Multichannel NMF

116[Ozerov&Févotte2010][Kameoka+2010][Sawada+2012]• Probability density function of microphone array observations
local Gaussian model
NMF model

log-likelihood

where

and

Assumption on mixing matrix• None [Ozerov+2010][Sawada+2012]
Applicable for underdetermined system• Invertible [Kameoka+2010][Kitamura+2015]
Specialized for determined system

Optimization• EM algorithm[Ozerov+2010][Kameoka+2010]• Auxiliary function approach[Sawada+2012][Kitamura+2015]

Multivariate extension of IS divergence• Different expression of: spatial property of source 𝑛𝑛• Multivariate extension of IS divergence(a.k.a. log-determinant divergence)

117Treating NMF bases as individual sources
Special case where
NMF bases:
Frequent spectral
patterns

allocated to eachseparation according
to the spatial property
of NMF bases

118Nonnegativity: from scalar to vector
ScalarVectorObservation
Nonnegative
representation
Low-rankmodel

nonnegativecomplex• Hermitian positive-semidefinite matrix

All eigenvalues are nonnegative• Spatial property of the k-th NMF basis at frequency i

119

120

Optimization problem of multichannel NMF• Objective function to be minimized
typically IS divergence-element wise:

withAuxiliary function design
Objective function[Sawada+2012, Higuchi+2014, Yoshii+2013]

withmatrix extension

Jensen's inequality applied
to reciprocal function

Auxiliary function

121

matrix extension

1st order Taylor expansionof logarithmic function

Multiplicative update rules
Multichannel IS-NMF

Algebraic Riccati equation

IS-NMF

122

123

Learned spatial property example
Basis-wise
The 10 bases seem to form 2 clusters,
each of which corresponds to each source.
source 1

source 2

Inter-channel phase difference of source 𝑘𝑘
becomes
where 𝐶𝐶 is speed of soundClustering NMF bases for sources• Modify the modeling of the spatial property•

with source-wise spatial property

Similar multiplicative updates derived

124

Learned spatial property example
Basis-wise

125

Source-wise

source 1

source 2

Inter-channel phase difference of source 𝑘𝑘
becomes
where 𝐶𝐶 is speed of sound3 music parts separation example

126• 3 sources and 2 microphones (underdetermined case)

The computational burden
was heavy: it took 838.30
seconds for separating 24second mixture.

4 examples in total can be found at
http://www.kecl.ntt.co.jp/icl/signal/sawada/demo/mchnmf/

Determined multichannel NMF

127[Kameoka+2010]• Special case of multichannel NMF where mixing matrix
is invertible• When, we can writeand sowhere

and• Earlier idea of Independent Low-Rank Analysis (ILRMA)BSS methods and optimization techniques
2000
Determined BSS

2010(Earlier idea of ILRMA)

ICA

Determined
MNMF

IVA

AuxIVA(MM)

FD-ICAILRMA(MM)TF maskestimation

Underdetermined BSS

MNMF
(EM)

MNMF
(MM)
MFHMM
(MM)

Monaural source
separation

NMF

IS-NMF

ComplexNMF

Factorial
HMM

Optimization techniques
Natural gradient

IS divergence minimization with MM

Separation matrix
optimization with IP

128

BSS methods and optimization techniques
2000
Determined BSS

2010(Earlier idea of ILRMA)

ICA

Determined
MNMF

IVA

AuxIVA(MM)

FD-ICAILRMA(MM)TF maskestimation

Underdetermined BSS

MNMF
(EM)

MNMF
(MM)
MFHMM
(MM)

Monaural source
separation

NMF

IS-NMF

ComplexNMF

Factorial
HMM

Optimization techniques
Natural gradient

IS divergence minimization with MM

Separation matrix
optimization with IP

129

Categorization of LGM-based BSS methods
Generative model
MethodAttias (2003)
Ozerov & Févotte (2010)
Duong et al. (2010)
Kameoka et al. (2010)
Yoshioka et al. (2011)
Ono et al. (2012)
Sawada et al. (2013)
Higuchi et al. (2014)
Kitamura et al. (2015)
López et al. (2015)
Adiloğlu & Vincent (2016)
Kounades-Bastian et al. (2016)130

Categorization of LGM-based BSS methods
Constraints on
MethodAttias (2003)
Ozerov & Févotte (2010)
Duong et al. (2010)
Kameoka et al. (2010)
Yoshioka et al. (2011)
Ono et al. (2012)
Sawada et al. (2013)
Higuchi et al. (2014)
Kitamura et al. (2015)
López et al. (2015)
Adiloğlu & Vincent (2016)
Kounades-Bastian et al. (2016)131

Tutorial structure
1. Introduction
1. Separation of audio/speech signals
2. Live demonstration
2. ICA and IVA
1. ICA: Independent Component Analysis
2. IVA: Independent Vector Analysis
3. NMF
1. NMF: Nonnegative Matrix Factorization
2. MNMF: Multichannel NMF
4. ILRMA
1. ILRMA: Independent Low-Rank Matrix Analysis

132

Tensor and sliced matrices
2.2. IVA
4.1. ILRMA

ChannelMixtureFrequency
Time

3.2. MNMF

133

Historical development(Over-)determined problem
2.1 Independent component analysis (ICA)
1994
1998

2.1 Frequency-domain ICA (FDICA)

Year

1999

2006
2009
2011

2016

3.1

Nonnegative matrix factorization (NMF)
Many applications of NMF
Generative models in NMF
Many extensions of NMF

2.2 Independent vector analysis (IVA)

2012
2013

Underdetermined problem

Many permutation solvers
were proposed for FDICA

2.2

3.1

Itakura–Saito NMF (ISNMF)

Auxiliary-function-based IVA (AuxIVA)
Time-varying Gaussian IVA

4.1

134

More precise spectral
model using NMF

Independent low-rank matrix analysis (ILRMA)

3.2

Multichannel NMF

More efficient calculation
for determined situation

Historical development(Over-)determined problem
2.1 Independent component analysis (ICA)
1994
1998

2.1 Frequency-domain ICA (FDICA)

Year

1999

2006
2009
2011

2016

3.1

Nonnegative matrix factorization (NMF)
Many applications of NMF
Generative models in NMF
Many extensions of NMF

2.2 Independent vector analysis (IVA)

2012
2013

Underdetermined problem

Many permutation solvers
were proposed for FDICA

2.2

3.1

Itakura–Saito NMF (ISNMF)

Auxiliary-function-based IVA (AuxIVA)
Time-varying Gaussian IVA

4.1

135

More precise spectral
model using NMF

Independent low-rank matrix analysis (ILRMA)

3.2

Multichannel NMF

More efficient calculation
for determined situation

Spatial and spectral models in BSS

136• Spatial model: assumption of mixing/demixing system• Spectral model: assumption for each source
SourcesignalMixingsystemSource Time-freq.
distribution structure

Observed
mixtureDemixing
systemEstimated
signalSpatial model- Instantaneous mixture in frequency domain- Liner time-invariant mixing system
Spectral model- Super Gaussianity of source dist.
- Sparse or low-rank time-freq. structure

137

IVA revisited: model• IVA extends ICA to multivariate probabilistic model[Hiroe+, 2006], [Kim+, 2006], [Kim+, 2007]

Frequency-domain ICA: frequency scalar random variables
IVA: frequency vector random variables
Mutually
independent

MixingmatrixObserved
Estimated
Demixing
vectorvectormatrixHigher-ordercorrelations
Permutation-problem-free estimation of

…
…

…

…

…

Spherical
distribution

Sourcevector

138

IVA revisited: compared to FDICA• Frequency-domain ICA (FDICA)

STFT

Update separation filter so that the estimatedsignals obey non-Gaussian distribution• IVA

Frequency

Demixing
systemSource obeys
non-Gaussian dist.
Estimation

Scalar r.v.s

Time

STFT

Update separation filter so that the estimatedsignals obey non-Gaussian multivariate distribution

Frequency

Demixing
systemTime
Frequency

Estimation

Mutually
independent

Time

Vector r.v.s

Observation

Non-GaussianCurrent
empirical dist. source dist.

Frequency

Mixture is close to Gaussian signal
because of central limit theorem
Observation

Time

Current
empirical dist.

Non-Gaussian
spherical
source dist.

Mutually
independent• Frequency vector spectral model (IVA)
Co-occurrence among frequency binsof each source

Extend vector model to
low-rank matrix model

Frequency

Extension of vector source model in IVA

Time

Vector-activated modelMore precise representation oftime-frequency structure

Frequency• NMF spectral model
Co-occurrence among time-frequency slotsof each source with a low-rank structure

139

Incrementation of frequency basesTime

Low-rank modelISNMF revisited: low-rank spectral model140• Itakura–Saito NMF (ISNMF) [Févotte+, 2009]

Minimization of the above is equivalent to a maximum likelihood(ML) estimation with a following generative model:
At each time-frequency slot, complex-valued component
VarianceIf

can be decomposed as, thenobeys

Parameters arealso decomposedISNMF revisited: low-rank spectral model141• Itakura–Saito NMF (ISNMF) [Févotte+, 2009]: Power spectrogram

Im

Small variance

Frequency bin

ReGrayscale shows
intensity of varianceTime frame
ReIm

Large varianceNOTE:
because the variancetime and frequency

is non-Gaussian
can fluctuate along

ILRMA: unified method of IVA and NMF

142• Independent low-rank matrix analysis (ILRMA) [Kitamura+, 2016]
Unification of
1. Estimation of demixing matrix
2. Low-rank approximation using
Observed
signalMaximize
independence

IVA(IVA spatial model)(ISNMF spectral model)

Estimated
signal Low-rank
modeling

ISNMF
ISNMF

Low-rank
spectral modelILRMA: unified method of IVA and NMF

143• Independent low-rank matrix analysis (ILRMA) [Kitamura+, 2016]
Motivation of ILRMA?
When sources are mixed…

Mix
Time

Time

Frequency
Frequency(2) Rank of spectrogram increases
Frequency
Frequency(1) Mixture dist. approaches to Gaussian

Time

Time

To separate the sources…
Maximize non-Gaussianity Restrict rank of separated signal(using NMF theory)(using ICA or IVA theory)
Low-rank assumption can avoid
the permutation problem

Cost function in ILRMA

144• Cost function in FDICA, IVA, or ILRMA

Source distribution(spectral model)

FDICA (Laplace)

Separable for frequency

e.g., [Sawada+, 2003]

IVA (spherical Laplace)[Hiroe+, 2006], [Kim+, 2006]

IVA (time-varying Gauss)[Ono+, 2012]

ILRMA (ISNMF model)[Kitamura+, 2016]

Frequency vector

Cost function in ILRMA

145• Cost function in ILRMA
Demixing matrix:

Estimated signal:

Cost function in
Cost function intime-varying Gaussian IVA
ISNMF(estimates demixing matrix) (estimates low-rank spectral structure)
Spatial model update: AuxIVA
Spectral model update: ISNMF

All the variables arealternatively updatedUpdate rule of parameters in ILRMA• Maximum-likelihood-based update rules

Alternatively update both models
Spatial model
Spectral model(demixing matrix)(NMF variables)

Update demixing filter [Ono+, 2011]

Update NMF parameters [Nakano+, 2010]

Update estimated signal

Update estimated variance: one-hot vector with one
Matlab codes for ILRMA:at th element
https://github.com/d-kitamura/ILRMA
See also (pseudo code)
http://d-kitamura.net/pdf/misc/AlgorithmsForIndependentLowRankMatrixAnalysis.pdf

146

Clustering NMF bases for sources• NMF bases should be automatically clustered into
each source
Introduce
partitioning function

Fixed number ofbases for each source

Adaptive number ofbases for each source

147

Update rule of parameters in ILRMA• Maximum-likelihood-based update rules with

Alternatively update both models
Spatial model
Spectral model(demixing matrix)(NMF variables)

Update demixing filter [Ono+, 2011]

Update NMF parameters [Nakano+, 2010]

Update estimated signal
Update estimated variance: one-hot vector with one
at th element
See also (pseudo code)
http://d-kitamura.net/pdf/misc/AlgorithmsForIndependentLowRankMatrixAnalysis.pdf

148

Historical development(Over-)determined problem
2.1 Independent component analysis (ICA)
1994
1998

2.1 Frequency-domain ICA (FDICA)

Year

1999

2006
2009
2011

2016

3.1

Nonnegative matrix factorization (NMF)
Many applications of NMF
Generative models in NMF
Many extensions of NMF

2.2 Independent vector analysis (IVA)

2012
2013

Underdetermined problem

Many permutation solvers
were proposed for FDICA

2.2

3.1

Itakura–Saito NMF (ISNMF)

Auxiliary-function-based IVA (AuxIVA)
Time-varying Gaussian IVA

4.1

149

More precise spectral
model using NMF

Independent low-rank matrix analysis (ILRMA)

3.2

Multichannel NMF

More efficient calculation
for determined situation

Multichannel NMF revisited: model• Multichannel NMF with full-rank covariance [Sawada+, 2013]
Multichannel
vector

Instantaneous spatial covarianceSpatial covariances in Spatial covarianceseach time-frequency slot of each source

Partitioning function
Basis matrix Activation matrix

Gains
Spectral patterns

Observed
Spatial model
Spectral model
multichannel signal Spatial property of each source Spectral patterns of all sources

150

ILRMA and multichannel NMF• Relationship b/w ILRMA and multichannel NMF?
Source distribution: same
Spatial model: different
ILRMA

Multichannel NMF

Instantaneous mixture
in frequency domain

Mixture of full-rank covariances(and power spectrograms)• Rank-1 spatial model [Duong+, 2010]
Equivalent to instantaneous mixture: steering vector(column vector of mixing matrix)

151

ILRMA and multichannel NMF• Multichannel NMF with rank-1 spatial model
Cost function of multichannel NMF

Substitute rank-1 spatial model

,

into

152

ILRMA and multichannel NMF

153• Multichannel NMF with rank-1 spatial model
Substituteinto the cost function of multichannel NMF

Transform the variables as

andCost in
ILRMASummary of ILRMA• From IVA side:
Introduce NMF spectral model (basis incrementation)• From multichannel NMF side:

LimitedSpatial model

Flexible

Introduce rank-1 spatial model (instantaneous mixture assumption)
Multichannel
NMF
Rank-1 spatial
model

IVA

NMF spectral
model

ILRMA

Limited Spectral model

Flexible154

Experiment• Conditions
Source signals
Analysis window
Shift length
Number of bases
Evaluation score

Music signals obtained from SiSEC2011
Two microphones and two sources (determined)
512-ms-long Hamming window
128 ms (1/4 shift)
30 per each source/60 for all sources
Improvement ot signal-to-distortion ratio (SDR)

Impulse response E2A(reverberation time: 300 ms)

155

156

Experiment• Two source case (ultimate nz tour)
Good

Poor

AuxIVA(Laplace)

Multichannel
NMF

ILRMA
(without Z)

ILRMA(with Z)

Experiment• Conditions
Source signals
Analysis window
Shift length
Number of bases
Evaluation score

Music signals obtained from SiSEC2011
Three microphones and three sources (determined)
512-ms-long Hamming window
128 ms (1/4 shift)
30 per each source/90 for all sources
Improvement ot signal-to-distortion ratio (SDR)

Impulse response E2A(reverberation time: 300 ms)

157

158

Experiment• Three source case (bearlin-roads, 14 s)
Good
ILRMA w/ Z, 60.7 s
ILRMA w/o Z, 15.1 s
IVA, 11.5 s

Multichannel NMF, 7647.3 s

Multichannel NMF
without Z
with Z

Other demo:
http://d-kitamura.net/en/demo_en.htm

Poor

159

Conclusion
Complex-valued Tensor• Multichannel audio mixture• Frequency domain via STFT

ICA, IVAILRMA• Independence• Super-Gaussian

NMF, MNMF• Spectral bases• Low-rank

Component → Vector → Low-rank Matrix

Auxiliary function-based optimization• Fast convergence• Problem dependent, but “cookbook” works

Advertisement: book chapters

160• MNMF and ILRMA will be published from Springer inMarch, 2018!
Audio Source Separation (Signals andCommunication Technology)
1st ed. 2018 by Ed. Shoji Makino
Ch. 5
General formulation of multichannel extensions of
NMF variants;
Hirokazu Kameoka, Hiroshi Sawada, and Takuya
Higuchi.
Ch. 6
Determined Blind Source Separation with
Independent Low-Rank Matrix analysis;
Daichi Kitamura, Nobutaka Ono, Hiroshi Sawada,
Hirokazu Kameoka, and Hiroshi Saruwatari.

161

Thank you very muchfor attending this tutorial !
Hiroshi Sawada

Nobutaka Ono

Hirokazu Kameoka

Daichi Kitamura

Supplement

Scope of this tutorial

Supplement p.2

Blind Audio Source Separation

Noover-determined

under-determined

ICA, IVA

Clustering(e.g. GMM)

ILRMA

Multi-channel NMF
Yes

Utilize training data

Multi-channel

Single-channel

NMF

DNN-based methods

Matlab codes for ILRMA: https://github.com/d-kitamura/ILRMA

Tutorial structure (detailed)
1. Introduction
1. Separation of audio/speech signals (pp. 3, 5-12)
2. Live demonstration (pp. 4)
2. ICA and IVA
1. ICA: Independent Component Analysis• Formulation and model (pp. 14-24)• Optimization methods (pp. 25-35)• Permutation and scaling problems (pp. 36-40)
2. IVA: Independent Vector Analysis• Formulation and source models (pp. 43-53)• Auxiliary function based optimization(pp. 54-62)• Experiments and faster algorithm in stereocase (pp. 63-74)

Supplement p.3

3. NMF
1. NMF: Nonnegative Matrix Factorization• Interpretation as generative model (pp. 77-85)• Optimization methods (pp. 86-107)• Expectation-Maximization algorithm (pp. 87, 88)• Auxiliary function approach (pp. 89-107)• Inequalities for auxiliary function design (pp. 99-105)

2. MNMF: Multichannel NMF• Formulation (pp. 110-119)• Auxiliary function based optimization(pp. 120-122)• Experiments (pp. 123-126)• Model categorization (pp. 127-131)
4. ILRMA
1. ILRMA: Independent Low-Rank Matrix Analysis• Developed from IVA (pp. 135-148)• Developed from MNMF (pp. 149-153)• Summary and experiments (pp. 154-158)

S ELECTED REFERENCES
Frequency-domain approach to convolutive BSS•[Smaragdis, 1998], [Parra and Spence, 2000], [Schobben and Sommen, 2002],[Anemüller and Kollmeier, 2000], [Asano et al., 2003], [Saruwatari et al., 2003],[Saruwatari et al., 2006], [Yoshioka et al., 2008], [Duong et al., 2010]• Permutation alignment [Sawada et al., 2004], [Sawada et al., 2007], [Sawada et al., 2011]• Scaling adjustment to microphone observations [Cardoso, 1998], [Murata et al., 2001],[Matsuoka and Nakashima, 2001], [Takatani et al., 2004], [Mori et al., 2006]
ICA: Independent Component Analysis•

Information-maximization approach [Bell and Sejnowski, 1995]• Maximum likelihood (ML) estimation [Cardoso, 1997]• Natural gradient [Amari et al., 1996], [Cichocki and Amari, 2002]• Equivariance property [Cardoso and Souloumiac, 1996]• FastICA [Hyvärinen et al., 2001]• Complex-valued ICA [Bingham and Hyvärinen, 2000], [Sawada et al., 2003]• Auxiliary function based ICA [Ono and Miyabe, 2010]
IVA: Independent Vector Analysis•

Multivariate p.d.f. [Hiroe, 2006], [Kim et al., 2006], [Kim et al., 2007]• FastIVA [Lee et al., 2006], [Lee et al., 2007]• HEAD problem [Yeredor, 2009]• Auxiliary function based IVA (AuxIVA) [Ono, 2011], [Ono, 2012b], [Ikeshita et al., 2017]• Time-varying Gaussian p.d.f. [Ono, 2012a], [Ono et al., 2012]• Supervised or model-based IVA [Ono et al., 2012], [Lopez et al., 2015],[Nesta and Koldovský, 2017]• Online IVA [Kim, 2010], [Taniguchi et al., 2014], [Sunohara et al., 2017]
NMF: Nonnegative Matrix Factorization•

Auxiliary function based optimization for Euclidean distance NMF and generalized
Kullback-Leibler divergence NMF [Lee and Seung, 1999], [Lee and Seung, 2001]• EM-based optimization for Itakura-Saito divergence NMF [Févotte et al., 2009]• Auxiliary function based optimization for Itakura-Saito divergence NMF [Kameoka et al., 2006],[Nakano et al., 2010], [Févotte and Idier, 2011]• Auxiliary function based optimization for β divergence NMF [Nakano et al., 2010],[Févotte and Idier, 2011]• Auxiliary function based optimization for sparse NMF [Kameoka et al., 2009]
Multi-channel NMF••

EM-based optimization [Ozerov and Févotte, 2010]
Auxiliary function based optimization [Sawada et al., 2012], [Sawada et al., 2013],[Higuchi and Kameoka, 2014]

ILRMA: Independent Low-Rank Matrix Analysis•

Earlier idea (determined multi-channel NMF) [Kameoka et al., 2010]
Multichannel NMF with rank-1 spatial model [Kitamura et al., 2015a]• ILRMA [Kitamura et al., 2016], [Kitamura et al., 2018]• Relaxation of Rank-1 spatial model [Kitamura et al., 2015b]• Maximization-equalization algorithm [Mitsui et al., 2017b]• Optimal window length [Kitamura et al., 2017]• Based on Student’s t-distribution [Mogami et al., 2017]• With sparse regularization [Mitsui et al., 2017a]• With spatial regularization [Mitsui et al., 2018]•

Other references related to auxiliary function based optimization•[Lange et al., 2000], [Hunter and Lange, 2000], [Hunter and Lange, 2004], [Ono et al., 2009],[Ono and Sagayama, 2010], [Févotte and Idier, 2011], [Yoshii et al., 2013],[Kameoka and Takamune, 2014], [Sun et al., 2017]
R EFERENCES[Amari et al., 1996] Amari, S., Cichocki, A., and Yang, H. H. (1996). A new learning algorithm for blind signal separation.
In Advances in Neural Information Processing Systems, volume 8, pages 757–763. The MIT Press.
[Anemüller and Kollmeier, 2000] Anemüller, J. and Kollmeier, B. (2000). Amplitude modulation decorrelation forconvolutive blind source separation. In Proc. ICA 2000, pages 215–220.
[Asano et al., 2003] Asano, F., Ikeda, S., Ogawa, M., Asoh, H., and Kitawaki, N. (2003). Combined approach of arrayprocessing and independent component analysis for blind separation of acoustic signals. IEEE Trans. Speech Audio
Processing, 11(3):204–215.
[Bell and Sejnowski, 1995] Bell, A. and Sejnowski, T. (1995). An information-maximization approach to blind separation andblind deconvolution. Neural Computation, 7(6):1129–1159.
[Bingham and Hyvärinen, 2000] Bingham, E. and Hyvärinen, A. (2000). A fast fixed-point algorithm for independent
component analysis of complex valued signals. International Journal of Neural Systems, 10(1):1–8.
[Cardoso, 1997] Cardoso, J.-F. (1997). Infomax and maximum likelihood for blind source separation. IEEE Signal Processing
Letters, 4(4):112–114.
[Cardoso, 1998] Cardoso, J.-F. (1998). Multidimensional independent component analysis. In Acoustics, Speech and Signal
Processing, 1998. Proceedings of the 1998 IEEE International Conference on, volume 4, pages 1941–1944. IEEE.
[Cardoso and Souloumiac, 1996] Cardoso, J.-F. and Souloumiac, A. (1996). Jacobi angles for simultaneous diagonalization.
SIAM Journal on Matrix Analysis and Applications, 17(1):161–164.
[Cichocki and Amari, 2002] Cichocki, A. and Amari, S. (2002). Adaptive Blind Signal and Image Processing. John Wiley &Sons.
[Duong et al., 2010] Duong, N., Vincent, E., and Gribonval, R. (2010). Under-determined reverberant audio source separation
using a full-rank spatial covariance model. IEEE Trans. Audio, Speech, and Language Processing, 18(7):1830–1840.
[Févotte et al., 2009] Févotte, C., Bertin, N., and Durrieu, J.-L. (2009). Nonnegative matrix factorization with the
Itakura-Saito divergence: With application to music analysis. Neural Computation, 21(3):793–830.
[Févotte and Idier, 2011] Févotte, C. and Idier, J. (2011). Algorithms for nonnegative matrix factorization with the
β-divergence. Neural computation, 23(9):2421–2456.
[Higuchi and Kameoka, 2014] Higuchi, T. and Kameoka, H. (2014). Joint audio source separation and dereverberation basedon multichannel factorial hidden markov model. In Machine Learning for Signal Processing (MLSP), 2014 IEEE
International Workshop on, pages 1–6. IEEE.
[Hiroe, 2006] Hiroe, A. (2006). Solution of permutation problem in frequency domain ICA using multivariate probability
density functions. In Proc. ICA 2006 (LNCS 3889), pages 601–608. Springer.
[Hunter and Lange, 2000] Hunter, D. R. and Lange, K. (2000). Quantile regression via an MM algorithm. Journal ofComputational and Graphical Statistics, 9(1):60–77.
[Hunter and Lange, 2004] Hunter, D. R. and Lange, K. (2004). A tutorial on mm algorithms. The American Statistician,58(1):30–37.
[Hyvärinen et al., 2001] Hyvärinen, A., Karhunen, J., and Oja, E. (2001). Independent Component Analysis. John Wiley &Sons.
[Ikeshita et al., 2017] Ikeshita, R., Kawaguchi, Y., Togami, M., Fujita, Y., and Nagamatsu, K. (2017). Independent vector
analysis with frequency range division and prior switching. In Signal Processing Conference (EUSIPCO), 2017 25th
European, pages 2329–2333. IEEE.
[Kameoka et al., 2006] Kameoka, H., Goto, M., and Sagayama, S. (2006). Selective amplifier of periodic and non-periodic
components in concurrent audio signals with spectral control envelopes. Proc. SIG Tech. Reports on Music and Computerof IPSJ, pages 77–84.
[Kameoka et al., 2009] Kameoka, H., Ono, N., Kashino, K., and Sagayama, S. (2009). Complex NMF: A new sparse
representation for acoustic signals. In Acoustics, Speech and Signal Processing, 2009. ICASSP 2009. IEEE International
Conference on, pages 3437–3440. IEEE.
[Kameoka and Takamune, 2014] Kameoka, H. and Takamune, N. (2014). Training restricted boltzmann machines with
auxiliary function approach. In Machine Learning for Signal Processing (MLSP), 2014 IEEE International Workshop on,
pages 1–6. IEEE.
[Kameoka et al., 2010] Kameoka, H., Yoshioka, T., Hamamura, M., Le Roux, J., and Kashino, K. (2010). Statistical model ofspeech signals based on composite autoregressive system with application to blind source separation. In International
Conference on Latent Variable Analysis and Signal Separation, pages 245–253. Springer.

[Kim, 2010] Kim, T. (2010). Real-time independent vector analysis for convolutive blind source separation. IEEE Trans.
Circuits and Systems I: Regular Papers, 57(7):1431–1438.
[Kim et al., 2007] Kim, T., Attias, H. T., Lee, S.-Y., and Lee, T.-W. (2007). Blind source separation exploiting higher-orderfrequency dependencies. IEEE Trans. Audio, Speech and Language Processing, 15(1):70–79.
[Kim et al., 2006] Kim, T., Eltoft, T., and Lee, T.-W. (2006). Independent vector analysis: An extension of ica to multivariatecomponents. In International Conference on Independent Component Analysis and Signal Separation, pages 165–172.
Springer.
[Kitamura et al., 2017] Kitamura, D., Ono, N., and Saruwatari, H. (2017). Experimental analysis of optimal window lengthfor independent low-rank matrix analysis. In Proc. EUSIPCO 2017, pages 1210–1214.
[Kitamura et al., 2015a] Kitamura, D., Ono, N., Sawada, H., Kameoka, H., and Saruwatari, H. (2015a). Efficient
multichannel nonnegative matrix factorization exploiting rank-1 spatial model. In Proc. ICASSP 2015, pages 276–280.
[Kitamura et al., 2015b] Kitamura, D., Ono, N., Sawada, H., Kameoka, H., and Saruwatari, H. (2015b). Relaxation of rank-1
spatial constraint in overdetermined blind source separation. In Proc. EUSIPCO 2015, pages 1271–1275.
[Kitamura et al., 2016] Kitamura, D., Ono, N., Sawada, H., Kameoka, H., and Saruwatari, H. (2016). Determined blindsource separation unifying independent vector analysis and nonnegative matrix factorization. IEEE/ACM Transactions on
Audio, Speech, and Language Processing, 24(9):1626–1641.
[Kitamura et al., 2018] Kitamura, D., Ono, N., Sawada, H., Kameoka, H., and Saruwatari, H. (2018). Determined blindsource separation with independent low-rank matrix analysis. In Makino, S., editor, Audio Source Separation. Springer.
[Lange et al., 2000] Lange, K., Hunter, D. R., and Yang, I. (2000). Optimization transfer using surrogate objective functions.
Journal of computational and graphical statistics, 9(1):1–20.
[Lee and Seung, 2001] Lee, D. and Seung, H. (2001). Algorithms for non-negative matrix factorization. In Advances in
Neural Information Processing Systems, volume 13, pages 556–562.
[Lee and Seung, 1999] Lee, D. D. and Seung, H. S. (1999). Learning the parts of objects with nonnegative matrix
factorization. Nature, 401:788–791.
[Lee et al., 2006] Lee, I., Kim, T., and Lee, T.-W. (2006). Complex FastIVA: A robust maximum likelihood approach of
MICA for convolutive BSS. In Proc. ICA 2006 (LNCS 3889), pages 625–632. Springer.
[Lee et al., 2007] Lee, I., Kim, T., and Lee, T.-W. (2007). Fast fixed-point independent vector analysis algorithms forconvolutive blind source separation. Signal Processing, 87(8):1859–1871.
[Lopez et al., 2015] Lopez, A. R., Ono, N., Remes, U., Palomaki, K., and Kurimo, M. (2015). Designing multichannel source
separation based on single-channel source separation. In Proc. ICASSP, pages 469–473.
[Matsuoka and Nakashima, 2001] Matsuoka, K. and Nakashima, S. (2001). Minimal distortion principle for blind source
separation. In Proc. ICA 2001, pages 722–727.
[Mitsui et al., 2017a] Mitsui, Y., Kitamura, D., Takamichi, S., Ono, N., and Saruwatari, H. (2017a). Blind source separation
based on independent low-rank matrix analysis with sparse regularization for time-series activity. In Proc. ICASSP 2017,
pages 21–25.
[Mitsui et al., 2017b] Mitsui, Y., Takamune, N., Kitamura, D., Saruwatari, H., Takahashi, Y., and Kondo, K. (2017b).
Independent low-rank matrix analysis based on parametric majorization-equalization algorithm. In Proc. CAMSAP 2017,
pages 98–102.
[Mitsui et al., 2018] Mitsui, Y., Takamune, N., Kitamura, D., Saruwatari, H., Takahashi, Y., and Kondo, K. (2018). Vectorwise
coordinate descent algorithm for spatially regularized independent low-rank matrix analysis. In Proc. ICASSP 2018.
[Mogami et al., 2017] Mogami, S., Kitamura, D., Mitsui, Y., Takamune, N., Saruwatari, H., and Ono, N. (2017). Independent
low-rank matrix analysis based on complex student’s t-distribution for blind audio source separation. In Proc. MLSP 2017.
[Mori et al., 2006] Mori, Y., Saruwatari, H., Takatani, T., Ukai, S., Shikano, K., Hiekata, T., Ikeda, Y., Hashimoto, H., andMorita, T. (2006). Blind separation of acoustic signals combining SIMO-model-based independent component analysis andbinary masking. EURASIP Journal on Applied Signal Processing, 2006:Article ID 34970, 17 pages.
[Murata et al., 2001] Murata, N., Ikeda, S., and Ziehe, A. (2001). An approach to blind source separation based on temporal
structure of speech signals. Neurocomputing, 41:1–24.
[Nakano et al., 2010] Nakano, M., Kameoka, H., Le Roux, J., Kitano, Y., Ono, N., and Sagayama, S. (2010).
Convergence-guaranteed multiplicative algorithms for nonnegative matrix factorization with β-divergence. In Machine
Learning for Signal Processing (MLSP), 2010 IEEE International Workshop on, pages 283–288. IEEE.
[Nesta and Koldovský, 2017] Nesta, F. and Koldovský, Z. (2017). Supervised independent vector analysis through pilot
dependent components. In Proc. ICASSP, pages 536–540.
[Ono, 2011] Ono, N. (2011). Stable and fast update rules for independent vector analysis based on auxiliary function
technique. In Applications of Signal Processing to Audio and Acoustics (WASPAA), 2011 IEEE Workshop on, pages
189–192. IEEE.
[Ono, 2012a] Ono, N. (2012a). Auxiliary-function-based independent vector analysis with power of vector-norm typeweighting functions. In Signal & Information Processing Association Annual Summit and Conference (APSIPA ASC), 2012
Asia-Pacific, pages 1–4.

[Ono, 2012b] Ono, N. (2012b). Fast stereo independent vector analysis and its implementation on mobile phone. In Proc.
International Workshop on Acoustic Echo and Noise Control (IWAENC). IEEE.
[Ono et al., 2009] Ono, N., Kohno, H., Ito, N., and Sagayama, S. (2009). Blind alignment of asynchronously recorded signals
for distributed microphone array. In Applications of Signal Processing to Audio and Acoustics (WASPAA), 2011 IEEE
Workshop on, pages 161–164. IEEE.
[Ono and Miyabe, 2010] Ono, N. and Miyabe, S. (2010). Auxiliary-function-based independent component analysis forsuper-Gaussian sources. In International Conference on Latent Variable Analysis and Signal Separation, pages 165–172.
Springer.
[Ono and Sagayama, 2010] Ono, N. and Sagayama, S. (2010). R-means localization: A simple iterative algorithm for source
localization based on time difference of arrival. In Proc. ICASSP, pages 2718–2721. IEEE.
[Ono et al., 2012] Ono, T., Ono, N., and Sagayama, S. (2012). User-guided independent vector analysis with source activity
tuning. In Acoustics, Speech and Signal Processing (ICASSP), 2012 IEEE International Conference on, pages 2417–2420.
IEEE.
[Ozerov and Févotte, 2010] Ozerov, A. and Févotte, C. (2010). Multichannel nonnegative matrix factorization in convolutive
mixtures for audio source separation. IEEE Trans. Audio, Speech and Language Processing, 18(3):550–563.
[Parra and Spence, 2000] Parra, L. and Spence, C. (2000). Convolutive blind separation of non-stationary sources. IEEE
Trans. Speech Audio Processing, 8(3):320–327.
[Saruwatari et al., 2006] Saruwatari, H., Kawamura, T., Nishikawa, T., Lee, A., and Shikano, K. (2006). Blind source
separation based on a fast-convergence algorithm combining ICA and beamforming. IEEE Trans. Audio, Speech andLanguage Processing, 14(2):666–678.
[Saruwatari et al., 2003] Saruwatari, H., Kurita, S., Takeda, K., Itakura, F., Nishikawa, T., and Shikano, K. (2003). Blindsource separation combining independent component analysis and beamforming. EURASIP Journal on Applied Signal
Processing, 2003(11):1135–1146.
[Sawada et al., 2007] Sawada, H., Araki, S., and Makino, S. (2007). Measuring dependence of bin-wise separated signals for
permutation alignment in frequency-domain BSS. In Proc. ISCAS 2007, pages 3247–3250.
[Sawada et al., 2011] Sawada, H., Araki, S., and Makino, S. (2011). Underdetermined convolutive blind source separation via
frequency bin-wise clustering and permutation alignment. IEEE Trans. Audio, Speech, and Language Processing,19(3):516–527.
[Sawada et al., 2012] Sawada, H., Kameoka, H., Araki, S., and Ueda, N. (2012). Efficient algorithms for multichannel
extensions of Itakura-Saito nonnegative matrix factorization. In Proc. ICASSP 2012, pages 261–264.
[Sawada et al., 2013] Sawada, H., Kameoka, H., Araki, S., and Ueda, N. (2013). Multichannel extensions of non-negative
matrix factorization with complex-valued data. IEEE Trans. Audio, Speech, and Language Processing, 21(5):971–982.
[Sawada et al., 2003] Sawada, H., Mukai, R., Araki, S., and Makino, S. (2003). Polar coordinate based nonlinear function forfrequency domain blind source separation. IEICE Trans. Fundamentals, E86-A(3):590–596.
[Sawada et al., 2004] Sawada, H., Mukai, R., Araki, S., and Makino, S. (2004). A robust and precise method for solving the
permutation problem of frequency-domain blind source separation. IEEE Trans. Speech Audio Processing, 12(5):530–538.
[Schobben and Sommen, 2002] Schobben, L. and Sommen, W. (2002). A frequency domain blind signal separation method
based on decorrelation. IEEE Trans. Signal Processing, 50(8):1855–1865.
[Smaragdis, 1998] Smaragdis, P. (1998). Blind separation of convolved mixtures in the frequency domain. Neurocomputing,
22:21–34.
[Sun et al., 2017] Sun, Y., Babu, P., and Palomar, D. P. (2017). Majorization-minimization algorithms in signal processing,communications, and machine learning. IEEE Transactions on Signal Processing, 65(3):794–816.
[Sunohara et al., 2017] Sunohara, M., Haruta, C., and Ono, N. (2017). Low-latency real-time blind source separation forhearing aids based on time-domain implementation of online independent vector analysis with truncation of non-causal
components. In Proc. ICASSP, pages 216–220.
[Takatani et al., 2004] Takatani, T., Nishikawa, T., Saruwatari, H., and Shikano, K. (2004). High-fidelity blind separation of
acoustic signals using SIMO-model-based independent component analysis. IEICE Trans. Fundamentals,
E87-A(8):2063–2072.
[Taniguchi et al., 2014] Taniguchi, T., Ono, N., Kawamura, A., and Sagayama, S. (2014). An auxiliary-function approach toonline independent vector analysis for real-time blind source separation. In Proc. HSCMA, pages 107–111.
[Yeredor, 2009] Yeredor, A. (2009). On hybrid exact-approximate joint diagonalization. In Computational Advances inMulti-Sensor Adaptive Processing (CAMSAP), 2009 3rd IEEE International Workshop on, pages 312–315. IEEE.
[Yoshii et al., 2013] Yoshii, K., Tomioka, R., Mochihashi, D., and Goto, M. (2013). Infinite positive semidefinite tensor
factorization for source separation of mixture signals. In International Conference on Machine Learning, pages 576–584.
[Yoshioka et al., 2008] Yoshioka, T., Nakatani, T., and Miyoshi, M. (2008). An integrated method for blind separation anddereverberation of convolutive audio mixtures. In Proc. EUSIPCO.