Algorithms for Independent Low-Rank Matrix Analysis
Daichi Kitamura1,*
1*

SOKENDAI (The Graduate University for Advanced Studies), Japand-kitamura@nii.ac.jp, http://d-kitamura.sakura.ne.jp/en/index en.html

ABSTRACTThis document summarizes an algorithm for independent low-rank matrix analysis, which was proposed as determined rank-1 multichannel nonnegative matrix factorization in the following published papers:
Daichi Kitamura, Nobutaka Ono, Hiroshi Sawada, Hirokazu Kameoka, and Hiroshi Saruwatari, “Efficient multichannel nonnegative
matrix factorization exploiting rank-1 spatial model,” Proceedings of IEEE International Conference on Acoustics, Speech and Signal
Processing (ICASSP 2015), pp. 276–280, April 2015.
Daichi Kitamura, Nobutaka Ono, Hiroshi Sawada, Hirokazu Kameoka, and Hiroshi Saruwatari, “Determined blind source separation
unifying independent vector analysis and nonnegative matrix factorization,” IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 24, no. 9, pp. 1626–1641, September, 2016 (open access).
Two detailed algorithms for implementation and some empirical knowledges for the use of them are described in this document.

1 Introduction
Nonnegative matrix factorization (NMF) [1]–[5] is a low-rank approximation (dimensionality reduction) with nonnegative
constraint, and many applications using NMF has been proposed, e.g., audio signal processing, text mining, image recognition,and bioinformatics. In particular, audio source separation is one of the most popular and successive research area of NMF [6]–[8]. For audio signals, NMF extracts some specific spectral patterns (bases) and their time-varying gains (activations) in theobserved signal. The source separation using NMF is achieved by clustering such decomposed nonnegative parts into eachsource, and this clustering may require some prior knowledge about the sources in the mixture (e.g., semi- and full-supervised
NMF methods [9, 10]).
As a blind source separation technique, multichannel NMF (MNMF) was proposed [11, 12], which simultaneously estimates the NMF variables (source model) and a spatial mixing matrix. The estimated bases and activations are attributed tosourcewise spatial features. These methods are generalized by Sawada et al. [13] exploiting more flexible spatial model called
full-rank spatial covariance [14]. In this method, similarly to the original MNMF, the estimated source model is clusteredbased on the estimated spatial model to achieve the separation. However, it is reported that these algorithms are sensitive tothe initial values of variables, and the computational cost is larger than those of conventional blind source separation methods.
For a determined or overdetermined blind source separation, independent component analysis (ICA) [15] and its extension
to the frequency domain (frequency-domain ICA: FDICA) [16]–[20] are popular and traditional techniques. These ICAbased methods estimate a demixing matrix (an inverse matrix of the mixing matrix) in the frequency domain. In particular,
independent vector analysis (IVA) [21]–[23] achieves better and stable source separation performance compared with FDICA,
and an auxiliary-function-based eﬃcient update rules for IVA were developed [24].
In [25, 26], a new algorithm for blind source separation was proposed, which is called determined rank-1 MNMF. Determined rank-1 MNMF has the same model as the conventional MNMF [11], but the demixing matrix is estimated insteadof the mixing matrix. This approach is similar to the traditional ICA-based blind source separation, and it is revealed that
determined rank-1 MNMF is a natural extension of IVA in terms of the flexibility of the source model (although the assumed
prior distributions in IVA and determined rank-1 MNMF are diﬀerent). Namely, the vector source model in IVA is extended
to the low-rank matrix source model using NMF decomposition in determined rank-1 MNMF. To clarify this issue, in thisdocument, we have renamed determined rank-1 MNMF as independent low-rank matrix analysis (ILRMA). The cost function
in ILRMA consists of those in IVA and simple NMF based on Itakura–Saito divergence. Therefore, all the variables including
demixing matrix and NMF variables can eﬃciently be optimized by an iterative projection proposed in [24] and auxiliaryfunction-based multiplicative update rules [4]. However, since these update rules derived in the paper [26] are described inscalar notations, it is diﬃcult to implement ILRMA in an eﬃcient way using matrix operations. In this document, the eﬃcientalgorithms for ILRMA using matrix notations are described in Sects. 2.2 and 2.3.

2 Algorithms for ILRMA
There are two models in ILRMA: ILRMA without partitioning function (hereafter referred to as ILRMA1) and ILRMA withpartitioning function (hereafter referred to as ILRMA2), where the partitioning function is used for determining the optimal
number of NMF bases for each source. Therefore, in ILRMA1, we must set the fixed number of bases for each source.
In ILRMA2, we only set the total number of bases, and the algorithm adaptively estimates the partition of all the bases by
optimizing the partitioning function. The detailed formulations of them can be found in [26].
2.1 Definitions

In this subsection, some definitions of variables are described as preliminaries. Note that the lowercase Italic characters
denote scalars, lowercase bold Italic characters denote column vectors, uppercase bold Italic characters denote matrices, anduppercase Sans-serif characters denote the third-order tensors. Also, T and H denote simple transpose and Hermitian transposeof vectors or matrices, respectively.
• M: number of channels (microphones) whose index is m
• N: number of sources whose index is n, where we assume M = N (overdetermined situation) in this document• I: number of frequency bins whose index is i
• J: number of time frames whose index is j
• L: number of bases for each source whose index is l
• K: total number of bases for all the sources whose index is k
• xi j = (xi j,1 xi j,2 · · · xi j,M )T : complex-valued observed (mixture) signals in time-frequency domain• X: complex-valued I × J × M tensor whose element is xi j,m• yi j = (yi j,1 yi j,2 · · · yi j,N )T : complex-valued estimated (separated) signals in time-frequency domain• Y: complex-valued I × J × N tensor whose element is yi j,n• Wi = (wi,1 wi,2 · · · wi,N )H : complex-valued N × M demixing matrix for the ith frequency bin, which leads yi j = Wi xi jH
xi j
• wi,n : complex-valued M × 1 vector in Wi called demixing filter for the nth source, where yi j,n = wi,n

• en : N × 1 vector with only the nth element equal to unity and the other elements equal to zeros• P: nonnegative I × J × N tensor that corresponds to the power spectrograms of estimated source signals• pi j,n : nonnegative element of P
• R: nonnegative I × J × N tensor that corresponds to the model spectrograms (time-frequency variances) for all sources(low-rank approximation of P using NMF decomposition)• Tn : nonnegative I × L basis matrix for the nth source, which is used in ILRMA1• Vn : nonnegative L × J activation matrix for the nth source, which is used in ILRMA1• til,n : nonnegative element of Tn .
• vl j,n : nonnegative element of Vn .
• T : nonnegative I × K basis matrix for all sources, which is used in ILRMA2• V : nonnegative K × J activation matrix for all source, which is used in ILRMA2• tik : nonnegative element of T
• vk j : nonnegative element of V
• Z = (z1 , z2 , · · · , zN )T : N × K matrix called partitioning function used in ILRMA2, where all the elements are in therange [0, 1]
2/7

• zn : K × 1 vector in Z
• znk : element of Z, where∑n znk = 1

• 1(size) : matrix of ones whose size is denoted as the superscript, e.g., 1(I×J)• ε: machine epsilon• ŷi j,n : complex-valued and scale-fitted M × 1 vector, which corresponds to the estimated nth source (spatial) image
Note that the third-order tensor with a subscript denotes the sliced matrix or the fiber vector in the original tensor [27]. For
example, since X is an I × J × M tensor, Xi:: denotes the J × M sliced matrix in X, X: j: denotes the I × M sliced matrix
in X, and X::m denotes the I × J sliced matrix in X. Also, Xi j: , Xi:m , and X: jm denote the M × 1, J × 1, and I × 1 fiber(column) vectors, respectively. In the algorithm description, ◦ and the quotient symbol for matrices denote the elementwise
multiplication and division, respectively. In addition, the absolute value | · | and the dotted exponent for matrices (e.g., X .2 )denote the elementwise absolute value and the elementwise exponent, respectively. For example, P::n = |Y::n |.2 . Moreover, themaximum operator max(·, ·) returns a matrix with the larger elements taken from two inputs in each entry.
2.2 Algorithm for ILRMA1

Since ILRMA1 does not utilize the partitioning function Z, the optimization is simply carried out by an alternative iterationof iterative-projection-based IVA update rule [24] and simple NMF update rules [3, 4]. Algorithm 1 shows a detailed process
flow in ILRMA1. For the implementation, we have to note the corresponding scalars, vectors, and matrices variables. For
H
example, when we update wi,n at the lines 13 and 14 in Algorithm 1, the nth row in Wi must be replaced to the new wi,nbefore it is used in the line 17.
2.3 Algorithm for ILRMA2

In ILRMA2, we only set the total number of bases. The K bases are shared (distributed) to each source by the partitioning
function Z. The optimization of Z is also carried out by an NMF-like multiplicative update rule. Algorithm 2 shows a detailed
process flow in ILRMA2.
2.4 Post process of ILRMA

The ICA-based source separation methods cannot determine scales of the estimated sources. Therefore, after the separation,a back-projection technique [28, 29] must be applied to the estimated signal yi j . This technique ideally restores the scalesof the estimated sources to their observed amplitudes in xi j , where the inverse matrix of the demixing matrix is utilized forthe restoration. Algorithm 3 shows the detailed process flow of the back-projection technique. The output signal ŷi j,n is ascale-fitted estimated signal of the nth source, where ŷi j,n is a multichannel estimated signal of the nth source (M × 1 vector),
which is often called the source image.

3 Empirical knowledges
In this section, we describe some empirical knowledges about the implementation and use of ILRMA.
• When the input multichannel signal xi j is overdetermined (M > N), principal component analysis should be applied inadvance for reducing its dimension so that M = N.
• When the input signal xi j includes zeros, the zero division might occur in NMF update rules. The zeros in xi j should
be replaced with the machine epsilons. Also, we must avoid to put zeros as the initial values of Tn , Vn , T , V , and Z
because the multiplicative update rules do not work for the zero elements.
• For the initial values of Wi , the identity matrix is preferable.
• When the number of time frames J is small, namely, the observed signal is too short, the matrix Wi Ui,n at the lines 13
in Algorithm 1 and 25 in Algorithm 2 tends to be a singular matrix, and its inversion cannot be calculated. In such acase, pseudo-inverse should be used, but the computational cost greatly increases in many languages (e.g., Matlab). We
can also avoid the inversion error by adding small values to the diagonal elements in the singular matrix.
• As a normalization coeﬃcient, we described the root mean powers of estimated sources at the lines 20 in Algorithm 1
and 33 in Algorithm 2. However, any kind of coeﬃcients can be used here. This normalization is processed to avoid adivergence of variables Wi or R to infinity. If the computational cost should be reduced, these normalization processes
can be omitted although the algorithms become unstable in some cases.
3/7

Algorithm 1: ILRMA1
input : observed multichannel complex-valued signals xi j
output: estimated complex-valued sources yi j
1 Initialize Wi with identity matrix or complex-valued random values and Tn and Vn with nonnegative random values;
2 Calculate yi j = Wi xi j for all i and j;// Initial estimated sources
3 Calculate P::n = |Y::n |.2 for all n;// Initial power spectrograms of estimated sources
4 Calculate R::n = Tn Vn for all n;// Initial model spectrograms
5 repeat6
for n = 1 to N do
7


[
]1 

(P ◦R.−2 )V T . 2 Tn ← max Tn ◦ ::nR.−1::nV T n, ε ;// Update of basis matrix

R::n = Tn Vn ;// New model spectrograms

::n

8


]. 1 [ T

T (P ◦R.−2 ) 2, ε;Vn ← max Vn ◦ n T T::nR.−1::n9n// New model spectrogramsfor i = 1 to I do
{ [
(
)]}T
(1×M)
Ui,n = 1J XHi:: Xi:: ◦ R.−1
; // Xi:: is J × M matrix, Ri:n is J × 1 vector, Ui,n is M × M matrix
i:n 1
wi,n ← (Wi Ui,n )−1 en ;// Update of demixing filter
1Hwi,n ← wi,n (wi,nUi,n wi,n )− 2 ;// Normalization of demixing filter
end

12
13
14
15

20

end
Calculate yi j = Wi xi j for all i and j;
Calculate P::n = |Y::n |.2 for all n;
for n = 1 todo
√ N∑
λn = IJ1 i, j pi j,n ;for i = 1 to I do
wi,n ← wi,n λ−1n ;end
P::n ← P::n λ−2n ;R::n ← R::n λ−2
n ;Tn ← Tn λ−2
n ;

21
22
23
24
25
26
27

// Update of activation matrix::n

R::n = Tn Vn ;

10
11

16
17
18
19

n

// New estimated sources// New power spectrograms of estimated sources// Normalization coefficient// Normalization of demixing filter// Normalization of separated power spectrogram// Normalization of model spectrogram// Normalization of basis matrix

end

28 until converge;

4/7

Algorithm 2: ILRMA2
input : observed multichannel complex-valued signals xi j
output: estimated complex-valued sources yi j
1 Initialize Wi with identity matrix or complex-valued random values, T and V with nonnegative random values, and Z with random

values in range [0, 1];().−1

2 Z ← Z ◦ 1(N×N) Z;
3 Calculate yi j = Wi xi j for all i and j;
4 Calculate P::n = |Y::n |.2 for all n;[()]
5 Calculate R::n = 1(I×1) znT ◦ T V for all n;
6 repeat7
for n = 1 to N do)1
( T(J×1) . 2{[T (P::n ◦R.−2::n )]◦V }1
;
8
an =
.−1
T
(J×1)
[(T R::n )◦V ]1
9
end
10
Z ← max (Z ◦ A, ε), where A = (a1 · · · aN )T ;
(
).−1
11
Z ← Z ◦ 1(N×N) Z ;
[(
)
]
12
Calculate R::n = 1(I×1) znT ◦ T V for all n;
13
for i = 1 to I do(
)1
T(N×1) . 2{[V (Pi:: ◦R.−2i:: )]◦Z }1
14
bi =
;
T
(N×1)
[(V R.−1i:: )◦Z ]1
15
end
16
T ← max (T ◦ B, ε), where B = (b1 · · · bI )T ;
[(
)
]
17
Calculate R::n = 1(I×1) znT ◦ T V for all n;
18
for j = 1 to J do)]}
 {[ T (. 12
T (N×1)
 T P: j: ◦R.−2
: j: ◦Z 1
[(
)
]

 ;
19
cj = 
T .−1
T (N×1)
T R: j: ◦Z

20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

1

∑
// Ensuring n znk = 1// Initial estimated sources// Initial power spectrograms of estimated sources// Initial model spectrograms// am is K × 1 vector// Update of partitioning function, A is N × K matrix∑// Ensuring n znk = 1// New model spectrogram// bi is K × 1 vector// Update of basis matrix, B is I × K matrix// New model spectrogram// c j is K × 1 vectorend
V ← max (V ◦ C, ε), where C = (c1 · · · c J );// Update of activation matrix, C is K × J matrix[()]
Calculate R::n = 1(I×1) znT ◦ T V for all n;// New model spectrogramfor n = 1 to N do
for i = 1 to I do
{ [
(
)]}T
(1×M)
Ui,n = 1J XHi:: Xi:: ◦ R.−1
; // Xi:: is J × M matrix, Ri:n is J × 1 vector, Ui,n is M × M matrix
i:n 1
wi,n ← (Wi Ui,n )−1 en ;// Update of demixing filter
1H
Ui,n wi,n )− 2 ;wi,n ← wi,n (wi,n// Normalization of demixing filter
end
end
Calculate yi j = Wi xi j for all i and j;// New estimated sources
Calculate P::n = |Y::n |.2 for all n;// New power spectrograms of estimated sources
for n = 1 to
do
√ N∑
1
λn = IJ i, j pi j,n ;// Normalization coefficientfor i = 1 to I do
wi,n ← wi,n λ−1n ;end
P::n ← P::n λ−2n ;R::n ← R::n λ−2n ;// Normalization of demixing filter// Normalization of separated power spectrogram// Normalization of model spectrogramend∑
Calculate tik ← tik n znk λ−2n for all i and k;

λ−2Calculate znk ← znk ∑ ′ zn′ λ−2 for all n and k;
n n k ′// Normalization of basis matrix// Normalization of partitioning functionn

42 until converge;

5/7

Algorithm 3: Back-projection techniqueinput : estimated multichannel complex-valued sources yi j
output: scale-fitted estimated source images ŷi j,n
1 for i = 1 to I do
2
for j = 1 to J do
3
for n = 1 to N do()
ŷi j,n = Wi−1 en ◦ yi j ;4// M × 1 vector

end

5

end

6
7 end

• The separation result depends on the initial values. The error bars sometimes become large especially for the mixture ofspeech signals (see the experimental results in the paper [26]). The dependence on initial values becomes strong when
the number of bases increases.
• Similarly to FDICA or IVA, ILRMA utilizes an assumption about the mixing system called the rank-1 spatial model(assumption of instantaneous mixture in the frequency domain). On the basis of this assumption, the reverberation timeshould be much shorter than the length of window function used in STFT. Also, the sources must spatially be stable(not moving) and must be assumed as point sources.

4 Conclusion
This document summarizes the algorithms for ILRMA. There are two models in ILRMA depending on the presence ofpartitioning function. I hope this document will help your happy source separation life.

5 Acknowledgements
I would like to thank Yoshiaki Bando and Yoshiki Mitsui for kindly implementing and checking the algorithms described inthis document.

References[1] D. D. Lee and H. S. Seung, “Learning the parts of objects by non-negative matrix factorization,” Nature, vol. 401, pp.
788–791, 1999.
[2] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” Proc. Advances in Neural Information
Processing Systems, vol. 13, pp. 556–562, 2000.
[3] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix factorization with the Itakura–Saito divergence. With application to music analysis,” Neural Computation, vol. 21, no. 3, pp. 793–830, 2009.
[4] M. Nakano, H. Kameoka, J. L. Roux, Y. Kitano, N. Ono, and S. Sagayama, “Convergence-guaranteed multiplicative algorithms for nonnegative matrix factorization with β-divergence,” Proc. IEEE International Workshop on Machine Learning
for Signal Processing, pp. 283–288, 2010.
[5] A. Cichocki, R. Zdunek, A. H. Phan, and S. Amari, “Nonnegative matrix and tensor factorizations: applications to
exploratory multi-way data analysis and blind source separation,” John Wiley & Sons, 2009.
[6] T. Virtanen, “Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria,” IEEE Trans. Audio, Speech, and Language Processing, vol. 15, no. 3, pp. 1066–1074, 2007.
[7] A. Ozerov, C. Févotte, and M. Charbit, “Factorial scaled hidden Markov model for polyphonic audio representation andsource separation,” Proc. IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, pp. 121–124,
2009.
[8] H. Kameoka, M. Nakano, K. Ochiai, Y. Imoto, K. Kashino, and S. Sagayama, “Constrained and regularized variantsof non-negative matrix factorization incorporating music-specific constraints,” Proc. IEEE International Conference on
Acoustics, Speech, and Signal Processing, pp. 5365–5368, 2012.
6/7

[9] P. Smaragdis, B. Raj, and M. Shashanka, “Supervised and semi-supervised separation of sounds from single-channel
mixtures,” Proc. International Conference on Independent Component Analysis and Signal Separation, pp. 414–421,
2007.
[10] D. Kitamura, H. Saruwatari, K. Yagi, K. Shikano, Y. Takahashi, and K. Kondo, “Music signal separation based onsupervised nonnegative matrix factorization with orthogonality and maximum-divergence penalties,” IEICE Trans. Fundamentals of Electronics, Communications and Computer Sciences, vol. E97-A, no. 5, pp. 1113–1118, 2014.
[11] A. Ozerov and C. Févotte, “Multichannel nonnegative matrix factorization in convolutive mixtures for audio source
separation,” IEEE Trans. Audio, Speech, and Language Processing, vol. 18, no. 3, pp. 550–563, 2010.
[12] H. Kameoka, T. Yoshioka, M. Hamamura, J. Le Roux, and K. Kashino, “Statistical model of speech signals based oncomposite autoregressive system with application to blind source separation,” Proc. International Conference on Latent
Variable Analysis and Signal Separation, pp. 245–253, 2010.
[13] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Multichannel extensions of non-negative matrix factorization withcomplex-valued data,” IEEE Trans. Audio, Speech, and Language Processing, vol. 21, no. 5, pp. 971–982, 2013.
[14] N. Q. K. Duong, E. Vincent, and R. Gribonval, “Under-determined reverberant audio source separation using a full-rank
spatial covariance model,” IEEE Trans. Audio, Speech, and Language Processing, vol. 18, no. 7, pp. 1830–1840, 2010.
[15] P. Comon, “Independent component analysis, a new concept?,” Signal Processing, vol. 36, no. 3, pp. 287–314, 1994.
[16] P. Smaragdis, “Blind separation of convolved mixtures in the frequency domain,” Neurocomputing, vol. 22, pp. 21–34,
1998.
[17] 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. Speech and Audio Processing, vol. 11, no. 2, pp.
109–116, 2003.
[18] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Convolutive blind source separation for more than two sources in
the frequency domain,” Proc. IEEE International Conference on Acoustics, Speech, and Signal Processing, pp. III-885–II
I-888, 2004.
[19] H. Buchner, R. Aichner, and W. Kellerman, “A generalization of blind source separation algorithms for convolutive
mixtures based on second order statistics,” IEEE Trans. Speech and Audio Processing, vol. 13, no. 1, pp. 120–134, 2005.
[20] H. Saruwatari, T. Kawamura, T. Nishikawa, A. Lee, and K. Shikano, “Blind source separation based on a fastconvergence algorithm combining ICA and beamforming,” IEEE Trans. Audio, Speech, and Language Processing, vol.
14, no. 2, pp. 666–678, 2006.
[21] T. Kim, T. Eltoft, and T.-W. Lee, “Independent vector analysis: an extension of ICA to multivariate components,” Proc.
International Conference on Independent Component Analysis and Blind Source Separation, pp. 165–172, 2006.
[22] A. Hiroe, “Solution of permutation problem in frequency domain ICA using multivariate probability density functions,”
Proc. International Conference on Independent Component Analysis and Blind Source Separation, pp. 601–608, 2006.
[23] T. Kim, H. T. Attias, S.-Y. Lee, and T.-W. Lee, “Blind source separation exploiting higher-order frequency dependencies,”
IEEE Trans. Audio, Speech, and Language Processing, vol. 15, no. 1, pp. 70–79, 2007.
[24] N. Ono, “Stable and fast update rules for independent vector analysis based on auxiliary function technique,” Proc. IEEE
Workshop on Applications of Signal Processing to Audio and Acoustics, pp. 189–192, 2011.
[25] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and H. Saruwatari, “Eﬃcient multichannel nonnegative matrix factorization exploiting rank-1 spatial model,” Proc. IEEE International Conference on Acoustics, Speech and Signal Processing,pp. 276–280, April, 2015.
[26] 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. Audio, Speech, and Language Processing, vol. 24, no. 9, pp. 1626–1641, September, 2016.
[27] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3, pp. 455–500,
2009.
[28] N. Murata, S. Ikeda, and A. Ziehe, “An approach to blind source separation based on temporal structure of speech
signals,” Neurocomputing, vol. 41, no. 1–4, pp. 1–24, 2001.
[29] K. Matsuoka and S. Nakashima, “Minimal distortion principle for blind source separation,” Proc. International Conference on Independent Component Analysis and Blind Signal Separation, pp. 722–727, 2001.

7/7