Supplementary Information
Error-Gated Hebbian Rule: A Local Learning Rule for Principal and Independent
Component Analysis
Takuya Isomura, Taro Toyoizumi

Supplementary Figures
Figure S1
Figure S2
Supplementary Movie
Movie S1
Supplementary Source Codes
scirep_fig2_eghr.cscirep_mathematics.h
scirep_functions.h
Supplementary Methods
S1. Lemma
S2. Fixed point of the EGHR-β
S3. Linear stability of the EGHR-β
S4. Conventional PCA algorithm

1Supplementary Figures

●
●
● ●
● ●
●● ●●
●
●
●
●●●●● ●
● ●●●●● ●●●●●●●●●●●●●●●● ●●● ●●●●●
●
●
●
●
●●● ●
● ●
●
●
●
●
●

●
●

2

106

107


Time [step]

106

Time [step]

107



8
14 20
Dimension

●

●

●●

●●●●
●
●
●
●
●
●● ●
● ●
●
●● ●●
●●
●
●
●
●
●
●●● ● ●
● ●●●●● ●●
●●●●
●●
●
●
●
●
●
● ●
●●●
● ●
●
●
●
●
●

2

8
14 20
Dimension

Figure S1. Dimension dependency of EGHR-β relaxation time with β = 0. The
relaxation time is the time needed to perform ICA, where we define that ICA isachieved when the ratio of first to second maximum absolute values for every row andcolumn of matrix K = WA is less than threshold eth = 0.1. The 2-to-20-dimensionalsources follow a unit Laplace or uniform distribution, and mixing matrix A is a squarematrix multiplying a distortion matrix and a random rotation matrix, where thedistortion matrix is defined such that the covariance of inputs becomes cov(xi, xi) = 1
and cov(xi, xj) = 0.4 for i ≠ j. The learning rate of η = 1 × 10–5 was used, and W wasstarted from an identity matrix. Simulations were conducted 10 times for eachdimension with different A. (A) Dimension dependency of the relaxation time with
Laplace sources. Red filled and open circles are the median and maximum relaxation
times, respectively. Yellow filled and open circles are the results of the original EHGR
[4] for comparison. Note that the Laplace prior distribution p0(si) ∝ exp(– 2 |si|) was
assumed to calculate g(u) and E(u). (B) Dimension dependency of the relaxation timewith uniform sources. Blue filled and open circles are the median and maximum
relaxation times, respectively. Green filled and open circles are the results of theoriginal EHGR [4]. Note that the prior distribution of p0(si) ∝ exp(–si4/4) was assumed
to calculate g(u) and E(u). Therefore, the EGHR-β as well as the original EGHR
reliably perform ICA with random mixing matrices and the up-to-20-dimensionalsources without being trapped in spurious solutions.
2105 106 107 108

Time [step]●●●●●●●●●●●●●●●●●●●●●●●●●●

●

●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●●●●●
●●●
●
●
●
●
●
●
●●
● ●
● ●●●
●●●

●
●

2−2

1

22

24 ∞

a
Figure S2. Robustness of EGHR-β with β = 0 to a choice of nonlinear function g(u).
Suppose, the two-dimensional sources follow p(si) ∝ exp(–b |si|a) (a > 0, b is definedsuch that the variance of si is one), and A = (cosπ/6, –sinπ/6; sinπ/6, cosπ/6) is a rotation
matrix. The same measure of relaxation time as Fig. S1 was plotted, but the threshold
was defined as eth = 0.01. The learning rate of η = 1 × 10–5 and the upper bound of
simulation time T = 1 × 108 were used, and W was started from an identity matrix. Red
and blue circles represent relaxation times with non-linear functions g(u) optimized fora Laplace distribution (thus a = 1; see the red arrow in the figure) and a uniform
distribution (a = ∞; blue arrow), respectively. Filled circles indicate that ICA wassuccessful with the non-linear function before the T-th step, while open circles indicatethat ICA was not achieved before the T-th step. Yellow and green circles are the resultsof the original EGHR [4] with non-linear functions g(u) optimized for a Laplace
distribution and a uniform distribution, respectively. Therefore, the EGHR-β as well as
the original EGHR robustly perform ICA for a range of a (except a = 2, where no ICA
solution exists) by using g(u) designed either for a = 1 or ∞.

3

Supplementary Movie
Supplementary Movie 1. Performance of EGHR-β with natural and noise images. See
the Figure 3 legend and Methods for details.

Supplementary source codes
Supplementary Source Codes 1-3. C-language source codes of EGHR-β that
demonstrate PCA and ICA. See the Figure 2 legend and Methods for details.

4Supplementary Methods
S1. Lemma
S1.1 When 〈●〉 = ∫ ● p0(s) ds and g(si) = –dlogp0(si)/dsi, 〈φ(s)g(s)〉 = 〈φ'(s)〉 holds forarbitrary function φ(s).
[Proof]Since g(s) = –p0'(s) / p0(s), we have〈φ(s)g(s)〉 = ∫ φ(s)g(s) p0(s) ds = –∫ φ(s) p0'(s) ds = ∫ φ'(s) p0(s) ds = 〈φ'(s)〉.
[Proof end]

S1.2 Suppose an M×M square matrix C. When s1, …, sM independently follow even
distributions p(s1), …, p(sM) with zero mean and unit variance, respectively, we have〈(sT C s) ssT〉 = C + CT + tr(C)I + Diag[κiCii]. Note that κi ≡ 〈si4〉 – 3.
[Proof]
Suppose G = 〈(sT C s) ssT〉. Since sT C s = ∑k ∑l Ckl sk sl, an element of G is Gij = 〈∑k∑l Ckl sk sl si sj〉. When i = j, since 〈si〉 = 〈si3〉 = 0, 〈si2〉 = 1, and 〈si4〉 = κi, + 3, we have
Gii = 〈∑k Ckk sk2 si2〉 = ∑k≠i Ckk + Cii〈si4〉 = 2Cii + κiCii + tr(C).
When i ≠ j, we have
Gij = 〈(Cij + Cji) si2 sj2〉 = Cij + Cji.
Therefore, we get G = C + CT + tr(C)I + Diag[κiCii].
[Proof end]

S1.3 Suppose an N×M rectangle matrix C, an M-dimensional vector s = (s1, …, sM)T,and an N-dimensional vector s! = (s1, …, sN)T. When s1, …, sN independently follow an
identical even distribution p0(si) with zero mean and unit variance, sN+1, …, sM
independently follow distributions with zero mean and unit variance, and E( s! ) – 〈E( s! )〉= –log p0( s! ) + 〈log p0( s! )〉 and g'(si) = ∂2E( s! )/∂si2 hold, we have 〈(E( s! )–〈E( s!
)〉)Diag[g'( s! )]CssT〉 = Ω○C, where ○ is Hadamard product, Ωii = cov(–log p0(si),g'(si)si2) for i = j, and Ωij = cov(–log p0(si), g'(si)) + cov(–log p0(sj), sj2)〈g'(si)〉Θ[j≤N] fori ≠ j. Note that Θ[j≤N] is 1 for j ≤ N and 0 otherwise.
5

[Proof]Since E( s! ) and g'(si) are even functions,M〈 (E( s! ) – 〈E( s! )〉) Diag[g'( s! )]CssT 〉ij = 〈 (E( s! ) – 〈E( s! )〉) g'(si) ∑ Cil sl s j 〉
l=1= 〈 (E( s! ) – 〈E( s! )〉) g'(si)sj2 〉Cij.
We define Ωij = 〈 (E( s! ) – 〈E( s! )〉) g'(si)sj2 〉. When i = j,NΩii = 〈 ∑ ( –log p0(sk)+〈log p0(sk)〉)g'(si)si2〉 = 〈(–log p0(si)+〈log p0(si)〉)g'(si)si2〉
k=1= cov(–log p0(si), g'(si)si2).
When i ≠ j,N

Ωij = 〈 ∑ ( –log p0(sk)+〈log p0(sk)〉)g'(si)sj2〉
k=1= 〈{–log p0(si)+〈log p0(si)〉 + (–log p0(sj)+〈log p0(sj)〉)Θ[j≤N]}g'(si)sj2〉= 〈(–log p0(si)+〈log p0(si)〉)g'(si)〉 + 〈(–log p0(sj)+〈log p0(sj)〉)sj2〉〈g'(si)〉Θ[j≤N]= cov(–log p0(si), g'(si)) + cov(–log p0(sj), sj2)〈g'(si)〉Θ[j≤N].
[Proof end]

S2. Fixed point of the EGHR-βSuppose that mixing matrix A ∈ ℝM×M consists of A = R Λ1/2 B without loss ofgenerality, where R and B ∈ ℝM×M are rotation matrices and Λ ∈ ℝM×M is a diagonal
matrix. Moreover, suppose s1, …, sM independently follow even distributions p(s1), …,
p(sM) with zero mean and unit variance. The cost function of the EGHR-β is given by
L≡α〈

1
β
γ
(E(u) – 〈E(u)〉)2 – E(u)〉 +〈(Eu – Ex)2〉 +tr(WTW),222where E(u), Ex and Eu are defined byE(u) ≡ –log p0(u),Eu ≡ –log N(u) + 〈log N(u)〉 =

1 2 1
1 T
1
|u| –
〈|u|2〉 =u u–tr(KTK),
2
2
2
2

6

(S1)

Ex ≡ –log N(x) + 〈log N(x)〉 =

1 2 1
1 T
1
|x| –
〈|x|2〉 =x x–tr(ATA).
2
2
2
2

(S2)

Note that α, β, and γ are constants that control fixed points and the linear stability, inwhich α + β = 1 holds without loss of generality. Partial derivative of L by W is givenby∂L= α〈(E(u) – E0) g(u) xT〉 + β〈(Eu – Ex) u xT〉 + γW
∂W

(S3)

and the first-order differential form of L is given by∂L
dKij
i=1 j=1 ∂ K ij
NM

dL = ∑ ∑NM= ∑ ∑{ α〈(E(u) – E0) g(ui) sj〉 + β〈(Eu – Ex) ui sj〉} dKij + γWij dWij
i=1 j=1= α〈(E(u) – E0) g(u)T du〉 + β〈(Eu – Ex) uT du〉 + γtr(WTdW).

(S4)

Note that g(ui) ≡ –dlogp0(ui)/dui and E0 ≡ 〈E(u)〉 + 1. Since 〈(sT C s) ssT〉 = C + CT +tr(C)I + Diag[κiCii] (where κi = 〈si4〉 – 3; see Supplementary Methods S1.2), we have〈(Eu – Ex) u xT〉 ==1〈( uTu – tr(KTK) – xTx + tr(ATA) ) u xT〉
2

1
K〈( sTKTKs – tr(KTK) – sTATAs + tr(ATA) ) s sT〉AT
2= K( KTK +1
1
Diag[κi(KTK)ii] – ATA –Diag[κi(ATA)ii] )AT2
2= K( Φ○(KTK) – Φ○(ATA) )AT,

(S5)

where ○ is Hadamard product and Φ ∈ ℝM×M is a symmetric matrix such that Φii ≡1+κi/2 for i = j and Φij ≡ 1 for i ≠ j. Thus, ∂L/∂W becomes∂L= α〈(E(u) – E0) g(u) sT〉AT + βK( Φ○(KTK) – Φ○(ATA) )AT + γKA–1A–TAT. (S6)∂W
Obviously, K = O is one of fixed points of Eq.(S6).
Case 1. Suppose B = I. When we use E(u) – 〈E(u)〉 = |u|2/2 – 〈|u|2〉/2, Eq.(S6) becomes∂L= αK( Φ○(KTK) – I )AT + βK( Φ○(KTK) – Φ○(ATA) )AT + γKA–1A–TAT∂W= K{ Φ○(KTK) – βΦ○(ATA) – αI + γ(ATA)–1 }AT7

= K{ Φ○(KTK) – βΦ○Λ – αI + γΛ–1 }AT.

(S7)

Since A is a regular matrix, from ∂L/∂W = O, using a vector 1 = (1, …, 1)T ∈ ℝM, anecessary and sufficient condition of a fixed point is given by
KKTK = K{(11T–Φ)○(KTK) + βΦ○Λ + αI – γΛ–1}.

(S8)

We define a diagonal matrix D ∈ ℝM×M such thatD ≡ (11T–Φ)○(KTK) + βΦ○Λ + αI – γΛ–1

⎛ D O ⎞
1⎟,
= Diag[–κi/2 (K K)ii + β(1+κi/2) Λii + α – γΛii ] = ⎜⎜⎝ O D2 ⎟⎠
T

–1

(S9)

where diagonal matrices D1 ∈ ℝN×N and D2 ∈ ℝ(M–N)×(M–N) are subsets of D. Moreover,K ∈ ℝN×M is represented as K = P(I, Q) without loss of generality, where P ∈ ℝN×N is aregular matrix and Q ∈ ℝN×(M–N) is any matrix. Thus, Eq.(S8) becomes⎛ D⎛I ⎞ T
1
P(I, Q) ⎜ T ⎟ P P(I, Q) = P(I, Q) ⎜
⎜ O
⎝Q ⎠
⎝

O ⎞
⎟,
D2 ⎟⎠

⟺ (I+QQT)PTP(I, Q) = (D1, QD2).

(S10)Each block becomes (I+QQT) PTP = D1 and (I+QQT) PTP Q = Q D2. By substituting theformer equation into the latter, we obtain D1 Q = Q D2. Since this is one of Sylvester
equations, it becomes (D1⊗I – I⊗D2)Vec[Q] = O, where ⊗ is the Kronecker productand Vec[Q] = (Q11, …, Q1(M–N), Q21, …, QN(M–N))T. If (D1⊗I – I⊗D2) is a regular matrix,
Q = O is a necessary and sufficient condition to satisfy the equation. Hereafter, suppose
Q = O. Then, a necessary and sufficient condition of fixed point is
K = (P, O),(S11)

where P is a regular matrix that holds
PTP = D1 = Diag[–κi/2 (PTP)ii + β(1+κi/2) Λii + α – γΛii–1],⟺ Φ1○(PTP) = βΦ1○Λ1 + αI – γΛ1–1,⟺ PTP = βΛ1 + Φ1†○(αI – γΛ1–1).

(S12)

Note that Φ1 ∈ ℝN×N is defined such that (Φ1)ij ≡ Φij for all 1 ≤ i, j ≤ N and Φ1† ∈ ℝN×N
8

is defined such that (Φ1†)ij ≡ (Φ1)ij–1 for all 1 ≤ i, j ≤ N. When (α, β, γ) = (0, 1, 0), P =JΛ11/2 with any rotation matrix J ∈ ℝN×N satisfies Eq.(S12). Thus, the EHGR obtain asubspace that N sources span in an equivalent way to the well-known Oja’s subspace
rule [11] except kurtosis term.
Case 2. Suppose (α, β, γ) = (1, 0, 0). Moreover, suppose s1, …, sN independently follow
an identical distribution p0(si), and sN+1, …, sM are independent of each other. Under thiscondition, K = (I, O) ∈ ℝN×M is an equilibrium point of the EHGR-β, where I is the
N×N identical matrix and O is the N×(M–N) zero matrix. This is confirmed bysubstituting K = (I, O) into Eq.(S6). Using an N-dimensional vector s! = (s1, …, sN)T,
E(u) and E0 become E( s! ) = –log p0( s! ) and E0 = 〈–log p0( s! )〉 + 1, respectively. Since〈φ(s)g( s! )〉 = 〈φ'(s)〉 for arbitrary function φ(s) from Supplementary Methods S1.1, we
have∂L= 〈(E( s! ) – 〈E( s! )〉 – 1) g( s! ) sT〉AT∂W= 〈g( s! ) sT + (E( s! ) – 〈E( s! )〉 – 1)(I, O)〉AT = O.

(S13)Therefore, K = (I, O) is a sufficient condition to be a fixed point of the EHGR-β. This
fixed point of the EGHR-β is equivalent to that of the original EGHR [4] while M canbe larger than N in this study.
Case 3. Suppose β and γ are small constants, s1, …, sN independently follow an
identical even distribution p0(si) with zero mean and unit variance, and sN+1, …, sM
independently follow distributions with zero mean and unit variance. Under thiscondition, as shown in the following, K = (I, O) + Δ = (I+Δ1, Δ2) ∈ ℝN×M is a sufficientcondition to be a fixed point of Eq.(S6), where I is the N×N identical matrix, O is theN×(M–N) zero matrix, Δ ≡ (Δ1, Δ2) ∈ ℝN×M is a matrix with small absolute value
elements, |Δij| << 1, and Δ1 ∈ ℝN×N and Δ2 ∈ ℝN×(M–N) are subsets of Δ. Using
N-dimensional vector s! = (s1, …, sN)T, g(u) becomes g(u) = g( s! +Δs) ≈ g( s! ) +Diag[g'( s! )]Δs. Since 〈(E( s! ) – 〈E( s! )〉 – 1) g( s! ) sT〉 = O and 〈φ(s)g( s! )〉 = 〈φ'(s)〉 forarbitrary function φ(s), the first term in the right side of Eq.(S6) becomes〈(E(u)–E0)g(u)sT〉 ≈ 〈{E( s! )–〈E( s! )〉+g( s! )TΔs–〈g( s! )TΔs〉–1}(g( s! )+Diag[g'( s! )]Δs)sT〉

9≈ 〈(g( s! )TΔs–〈g( s! )TΔs〉) g( s! ) sT〉 + 〈(E( s! )–〈E( s! )〉–1) Diag[g'( s! )]Δs sT〉= 〈(Diag[g'( s! )]Δs+(I,O)ΔTg( s! )) sT〉 + 〈g( s! )TΔs–〈g( s! )TΔs〉〉(I,O)+ 〈(E( s! )–〈E( s! )〉–1) Diag[g'( s! )]Δs sT〉= 〈(I,O)ΔTg( s! ) sT〉 + 〈(E( s! )–〈E( s! )〉) Diag[g'( s! )]Δs sT〉= (I,O)ΔT(I,O) + Ω○Δ= (Δ1T+Ω1○Δ1, Ω2○Δ2).

(S14)

Note that from Supplementary Methods S1.3, Ω ≡ (Ω1, Ω2) ∈ ℝN×M is defined such thatΩii ≡ cov(–log p0(si), g'(si)si2) for i = j and Ωij ≡ cov(–log p0(si), g'(si)) + cov(–log p0(sj),sj )〈g'(si)〉Θ[j≤N] for i ≠ j (Θ[j≤N] is 1 for j ≤ N and 0 otherwise). Note that Ω1 ∈ ℝ
2

N×N

and Ω2 ∈ ℝN×(M–N) are subsets of Ω. We define H ≡ (H1, H2; H3, H4) = ATA ∈ ℝM×M.
Since Δ, β, and γ are small, the second and third terms in the right side of Eq.(S6) isapproximated as β(I,O)( Φ○((I,O)T(I,O)) – Φ○(ATA) )AT + γ(I,O)(ATA)–1AT =(I,O)(βΦ○(I–H) + γH–1)AT. Therefore, Eq.(S6) becomes∂L≈ (Δ1T+Ω1○Δ1, Ω2○Δ2)AT – (I,O)(βΦ○(H–I) – γH–1)AT∂W= [(Δ1T+Ω1○Δ1, Ω2○Δ2)⎛ βΦ ! (H − I ) − γ (H −1 )⎞β H 2 − γ (H −1 )2
1
1
1⎟ ] AT– (I,O) ⎜−1−1⎜⎝
β H 3 − γ (H )3βΦ 2 ! (H 4 − I ) − γ (H )4 ⎟⎠= [(Δ1T+Ω1○Δ1, Ω2○Δ2) – (βΦ1○(H1 – I) – γ(H–1)1, βH2 – γ(H–1)2)] AT.

(S15)

To be an equilibrium point, Δ needs to satisfy
Δ1T+Ω1○Δ1 = βΦ1○(H1 – I) – γ(H–1)1andΩ2○Δ2 = βH2 – γ(H–1)2.

(S16)

Since the right side of the first equation is a symmetric matrix, if and only if Δ1 is asymmetric matrix, a solution exists. Suppose Δ1 is a symmetric matrix. Therefore, using
a vector 1 = (1, …, 1)T ∈ ℝN, we obtain
Δ1 = (11T+Ω1)†○{βΦ1○(H1 – I) – γ(H–1)1},Δ2 = Ω2†○{βH2 – γ(H–1)2}.

(S17)

Therefore, K = (I+Δ1, Δ2) is a sufficient condition of a fixed point. Since |Δij| << 1, this
fixed point is sufficiently close to s! = (s1, …, sN)T, thus, the EGHR-β can perform
PCA and ICA simultaneously. Specifically, if B = I, Δ becomes a diagonal matrix and10its (i, i) element is Δii = (1+Ωii)–1 {β Φii (Λii–1) – γ Λii–1}; thus, an exact solution of PCA
and ICA exists.
To see how the shape of the source distribution affects on the fixed point, suppose
p0(si) (1 ≤ i ≤ N) follow p0(si) = exp(–b|si|a)/Z, where a > 0 is a positive constant, b > 0 is
tuned in a manner to hold 〈si2〉 = 1, and Z is a partition function such that Z =∫exp(–b|si|a)dsi. Note that sN+1, …, sM follow distributions with zero mean and unit
variance. Since –log p0(si) = b|si|a + const., g(si) = –dlogp0(si)/dsi = ba|si|a–1sgn(si), andg'(si) = –dg(si)/dsi = ba(a–1)|si|a–2 + ba|si|a–1δ(si), we have –log p0(si) = g(si)si/a + const.,
g'(si)si = (a–1)g(si) and –log p0(si) g'(si) = (a–1)/a g(si)2 + const. g'(si). Thus, Ωij
becomes
For i = j,
Ωii = cov(–log p0(si), g'(si)si2) = 〈(–log p0(si)+〈log p0(si)〉) g'(si)si2〉= 〈(g(si)si/a – 〈g(si)si/a〉) (a–1)g(si)si〉= (a–1)/a 〈(g'(si)si+g(si)) si〉 + (a–1)/a 〈g(si)si – 〈g(si)si〉〉= (a–1)/a 〈(a–1)g(si)si + g(si)si〉 = a–1,
For i ≠ j ≤ N,
Ωij = cov(–log p0(si), g'(si)) + cov(–log p0(sj), sj2)〈g'(si)〉= 〈(–log p0(si)+〈log p0(si)〉) g'(si)〉 + 〈(–log p0(sj)+〈log p0(sj)〉) sj2〉 〈g'(si)〉= 〈(a–1)/a g(si)2 – 〈g(si)si/a〉 g'(si)〉 + 〈(g(sj)sj/a – 〈g(sj)sj/a〉) sj2〉 〈g'(si)〉= (a–2)/a 〈g'(si)〉 + 2/a 〈g'(si)〉 = 〈g'(si)〉,
For N+1 ≤ j ≤ M,
Ωij = cov(–log p0(si), g'(si)) = (a–2)/a 〈g'(si)〉.

11

(S18)

S3. Linear stability of the EGHR-β
The second derivative of L is given by

∂ dL
dKkl
k=1 l=1 ∂ K kl
NM

d2L = ∑ ∑= α 〈 (g(u)T du – 〈g(u)T du〉) g(u)T du + (E(u) – E0) duT Diag[g'(u)] du 〉
+ β 〈 (uT du – 〈uT du〉) uT du + (Eu – Ex) duT du 〉 + γ tr(dWT dW)= α 〈(sTdKTg(u) – 〈sTdKTg(u)〉) sTdKTg(u) + (E(u) – E0) sT dKT Diag[g'(u)] dK s〉+ β 〈(sTdKTKs – 〈sTdKTKs〉) sTdKTKs + (Eu–Ex) sTdKTdKs〉 + γ tr(A–TdKTdKA–1)= α tr[ 〈(sTdKTg(u) – 〈sTdKTg(u)〉)g(u)sT + (E(u) – E0) Diag[g'(u)] dK ssT〉dKT ]+ β tr[ 〈(sT dKT K s – 〈sT dKT K s〉)ssT〉dKTK + 〈(Eu – Ex)ssT〉dKTdK ]+ γ tr[A–1A–TdKTdK].

(S19)

From Supplementary Methods S1.2, 〈 (sT dKT K s – 〈sT dKT K s〉)ssT 〉 = dKTK + KTdK +Diag[κi(dKTK)ii] = Φ○(dKTK + KTdK) (where Φii = 1+κi/2 for i = j and Φij = 1 for i ≠ j).
Moreover, by Eq.(S5), 〈(Eu – Ex) ssT〉 = KTK + Diag[κi(KTK)ii]/2 – ATA –Diag[κi(ATA)ii]/2 = Φ○(KTK – ATA). Thus, the second term of Eq.(S19) becomestr[ 〈(sT dKT K s – 〈sT dKT K s〉)ssT〉dKTK + 〈(Eu – Ex)ssT〉dKTdK ]= tr[ {Φ○(dKTK + KTdK)} dKTK + {Φ○(KTK – ATA)} dKTdK ].

(S20)

Therefore, d2L becomes
d2L = α tr[ 〈(sTdKTg(u) – 〈sTdKTg(u)〉)g(u)sT + (E(u)–E0) Diag[g'(u)] dK ssT〉dKT ]+ β tr[{Φ○(dKTK+KTdK)} dKTK] + tr[{βΦ○(KTK–ATA) + γA–1A–T} dKTdK].

(S21)

Case 1. Suppose B = I. When we use E(u) – 〈E(u)〉 = |u|2/2 – 〈|u|2〉/2, K = (P, O) gives
an equilibrium point of the cost function wherever P satisfies Eq.(S12). The first term ofEq.(S21) becomestr[〈(sTdKTg(u) – 〈sTdKTg(u)〉)g(u)sT + (E(u)–E0) Diag[g'(u)] dK ssT〉dKT]= tr[〈(sTdKTu – 〈sTdKTu〉)usT + (uTu/2 – tr(KTK)/2 – 1) dK ssT〉dKT]= tr[K〈(sT dKTK s – tr(dKTK))ssT〉dKT + dK〈(sTKTKs/2 – tr(KTK)/2 – 1)ssT〉dKT]= tr[K(dKTK + KTdK + Diag[κi(dKTK)ii])dKT + dK(KTK + Diag[κi(KTK)ii]/2 – I)dKT]= tr[ Φ○(dKTK + KTdK) dKTK + (Φ○(KTK) – I) dKTdK ].

12

(S22)

Thus, by supposing K = (P, O) and dK = (dK1, dK2), d2L becomes
d2L = tr[(Φ○(dKTK))dKTK] + tr[(Φ○(KTdK))dKTK]+ tr[(Φ○(KTK – βATA) – αI + γA–1A–T)dKTdK]= tr[(Φ1○(dK1TP))dK1TP] + tr[(Φ1○(PTdK1),PTdK2)(PTdK1,PTdK2)T]⎛ Φ ! (PT P − βΛ ) − α I + γΛ −1
⎞
O
1
1
1
⎟ dKTdK]+ tr[ ⎜
−1⎜⎝O− βΦ 2 ! Λ 2 − α I + γΛ 2 ⎟⎠= tr[(Φ1○(dK1TP + PTdK1))dK1TP] + tr[(Φ1○(PTP–βΛ1) – αI + γΛ1–1)dK1TdK1]+ tr[PTdK2dK2TP] + tr[(–βΦ2○Λ2 – αI + γΛ2–1)dK2TdK2].

(S23)

Since Φ1○(PTP) = βΦ1○Λ1 + αI – γΛ1–1 (see Eq.(S12)) and Φ1○(dK1TP + PTdK1) =dK1TP + PTdK1 + Diag[κi(dK1TP)ii], d2L becomesN

d2L = tr[(dK1TP + PTdK1)dK1TP] + ∑ κi(dK1TP)ii2
i=1+ tr[PTdK2dK2TP] + tr[(–βΦ2○Λ2 – αI + γΛ2–1)dK2TdK2]
N= tr[(dX1T + dX1)dX1T] + ∑ κi(dX1T)ii2
i=1+ tr[dX2dX2T] + tr[(–βΦ2○Λ2 – αI + γΛ2–1)dX2T(P–1P–T)dX2]
=

N
1 N
2
(dX
+
dX
)
+
2
ij
ji
∑ Φii(dX1)ii2∑∑
2 i=1 j≠i≤Ni=1+ tr[dX2dX2T] – tr[(PTP)–1dX2(βΦ2○Λ2 – αI + γΛ2–1)dX2T].

(S24)

Note that dX = (dX1, dX2) = PTdK = PT(dK1, dK2) is supposed. Because coefficients of
dX1 are positive, the linear stability depends on the sign of coefficients of dX2. Since
PTP = βΛ1 + Φ1†○(αI – γΛ1–1) and (βΦ2○Λ2 – αI + γΛ2–1) are diagonal matrices, wehavetr[dX2dX2T] – tr[(PTP)–1dX2(βΦ2○Λ2 – αI + γΛ2–1)dX2T]
N

M

= ∑ ∑

i=1 j= N +1
NNM

dXij2 – ∑ ∑ ( (PTP)–1)ii dXij2 (βΦ2○Λ2 – αI + γΛ2–1)jj
i=1 j= N +1M

= ∑ ∑ ( (PTP)–1)ii dXij2 { (βΛ1+Φ1†○(αI–γΛ1–1))ii – (βΦ2○Λ2–αI+γΛ2–1)jj }i=1 j= N +1

13

N

M

= ∑ ∑ ( (PTP)–1)ii dXij2 { βΛii + Φii–1 (α–γΛii–1) – (βΦjjΛjj – α + γΛjj–1) }i=1 j= N +1
NM

= ∑ ∑ ( (PTP)–1)ii dXij2 { β(Λii–ΦjjΛjj) + α(Φii–1–1) – γ(Φii–1Λii–1–Λjj–1) }(S25)

i=1 j= N +1

Therefore, a necessary and sufficient condition for linear stability isβ(Λii–ΦjjΛjj) + α(Φii–1–1) – γ(Φii–1Λii–1–Λjj–1) > 0 for 1≤i≤N, N+1≤j≤M.

(S26)

When (α, β, γ) = (0, 1, 0) and sN+1, …, sM follow a unit Gaussian distribution, thecondition for linear stability is Λii > Λjj for 1 ≤ i ≤ N and N+1 ≤ j ≤ M. Thus, when u
represents a space that the first to Nth components span, the state is stable, while whenu involves other components, the state is unstable. Accordingly, the EGHR-β can
extract principal components. However, when sN+1, …, sM follow non-Gaussian
distributions, the linear stability condition depends on kurtosis of the distributions.
Case 2. Suppose (α, β, γ) = (1, 0, 0), s1, …, sN independently follow an even identical
distribution p0(si) with zero mean and unit variance, and sN+1, …, sM independently
follow distributions with zero mean and unit variance. Under this condition, K = (I, O)
is an equilibrium point. When K = (I, O), d2L becomes
d2L = tr[〈(sTdKTg( s! )–〈sTdKTg( s! )〉)g( s! )sT + (E( s! )–〈E( s! )〉–1) Diag[g'( s! )] dK ssT〉dKT]= tr[〈((I,O)dKTg( s! )+Diag[g'( s! )]dKs)sT + (E( s! )–〈E( s! )〉–1) Diag[g'( s! )] dK ssT〉dKT]= tr[ (I,O)dKT(I,O)dKT + 〈(E( s! )–〈E( s! )〉) Diag[g'( s! )] dK ssT〉dKT ]= tr[ dK1dK1 + (Ω○dK) dKT ]
N

N

N

M

= ∑ ( 1+Ωii)dKii + ∑ ∑ ( dKijdKji + ΩijdKij2) + ∑ ∑ Ωij dKij22

i=1

i=1 i≠ j≤NNN= ∑ ( 1+Ωii)dKii2 +i=1

i=1 j= N +1

NM1
( ΩijdKij2+2dKijdKji+ΩjidKji2) + ∑ ∑ Ωij dKij2.
∑∑
2 i=1 i≠ j≤N
i=1 j= N +1

(S27)

Therefore, necessary and sufficient conditions to be linearly stable are given by
1 + Ωii > 0

for 1 ≤ i = j ≤ N,

Ωij Ωji > 1

for 1 ≤ i ≠ j ≤ N,

Ωij > 0for 1 ≤ i ≤ N, N+1 ≤ j ≤ M.

(S28)

To see how the shape of the source distribution affects on the linear stability, again

14

suppose p0(si) (1 ≤ i ≤ N) to follow p0(si) = exp(–b|si|a)/Z and solve Ωij as shown in
Eq.(S18). Then, whenever a ≥ 1, 1 + Ωii = a > 0 and Ωij Ωji = 〈g'(si)〉2 > 1 hold. Whereas,Ωij = (a–2)/a 〈g'(si)〉 > 0 holds only when a > 2. Therefore, a > 2 is a necessary andsufficient condition to be linearly stable. Accordingly, when s1, …, sN follow a
sub-Gaussian distribution (a > 2) and sN+1, …, sM follow Gaussian or super-Gaussian
distributions, s1, …, sN are chosen as outputs. Whereas, when s1, …, sN follow asuper-Gaussian distribution (a < 2) and sN+1, …, sM follow Gaussian or sub-Gaussian
distributions, some of s1, …, sN may not be chosen while some of sN+1, …, sM may bechosen as outputs.
Case 3. Suppose that mixing matrix A ∈ ℝN×N consists of A = R Λ1/2, where R ∈ ℝN×N
is a rotation matrix and Λ ∈ ℝN×N is a diagonal matrix. Then, H becomes H = ATA ≈ Λ.
Moreover, suppose |β|, |γ| << 1 and α = 1 – β ≈ 1. From the fixed point analysis, when B= I, a fixed point is represented as u = (I+Δ1, O) s, where Δ1 ∈ ℝN×N is a diagonal
matrix defined by Δ1 = (11T+Ω1)†○{βΦ1○(Λ1 – I) – γΛ1–1}, while Δ2 = O. The first termof Eq.(S21) becomes〈(sT dKT g(u) – 〈sT dKT g(u)〉)g(u)sT + (E(u)–E0) Diag[g'(u)] dK ssT〉ij= 〈(∑k,l g(uk) dKkl sl – 〈∑k,l g(uk) dKkl sl〉) g(ui) sj + (E(u)–E0) g'(ui) ∑l dKil sl sj〉= 〈(∑k,l g(uk) dKkl sl – 〈∑k,l g(uk) dKkl sl〉) g(ui) sj + (E(u)–E0) g'(ui) dKij sj2〉.

(S29)

To further calculate Eq.(S29), again suppose p0(si) (1 ≤ i ≤ N) follow p0(si) =exp(–b|si|a)/Z. Thus, g(ui) becomes g(ui) = (1+Δii)a–1g(si) ≈ (1+(a–1)Δii)g(si) = Dii g(si),where D is diagonal matrix with Dii ≡ 1 + (a–1)Δii. Using this, we have〈∑k,l (g(uk) dKkl sl) g(ui) sj〉 ≈ Dii 〈∑k,l (Dkk g(sk) dKkl sl) g(si) sj〉= Dii 〈{∑l (Dii g'(si) dKil sl) + ∑k (Dkk g(sk) dKki)} sj + δij ∑k,l (Dkk g(sk) dKkl sl)〉
= Dii { Dii 〈g'(si)sj2〉 dKij + Djj dKji Θ[j≤N] + δij ∑k Dkk dKkk }= Dii2 〈g'(si)sj2〉 dKij + Dii Djj dKji Θ[j≤N] + Dii δij tr(DdK1).

(S30)

Moreover, by defining e(ui) ≡ –log p0(ui) + 〈log p0(ui)〉, we have〈(E(u)–E0) g'(ui) dKij sj2〉 = 〈(e(ui) + (1–δij)e(uj)Θ[j≤N] – 1) g'(ui) sj2〉 dKij
≈ 〈{(1+aΔii)e(si) + (1–δij)(1+aΔjj)e(sj)Θ[j≤N] – 1} (1+(a–2)Δii)g'(si) sj2〉 dKij
≈ 〈{(Dii+Δii)e(si) + (1–δij)(Dii+Δii–aΔii+aΔjj)e(sj)Θ[j≤N] – 1}(Dii–Δii)g'(si)sj2〉 dKij15

= {(Dii+Δii)(Dii–Δii)Ωij + (Dii–Δii)〈{(1–δij)a(Δjj–Δii)e(sj)Θ[j≤N] – 1}g'(si)sj2〉} dKij
≈ {Dii2 Ωij + 〈{(1–δij)a(Δjj–Δii)e(sj)Θ[j≤N] – (Dii–Δii)}g'(si)sj2〉} dKij
= {Dii2 Ωij + a(Δjj–Δii)〈g'(si)〉〈e(sj)sj2〉Θ[j≤N] – (Dii–Δii)〈g'(si)sj2〉} dKij
= {Dii2 Ωij + 2(Δjj–Δii)〈g'(si)〉Θ[j≤N] – (Dii–Δii)〈g'(si)sj2〉} dKij.

(S31)

In the last line, 〈e(sj)sj2〉 = 〈(g(sj)sj – 〈g(sj)sj〉)sj2〉/a = 2/a is used. Thus, by defining dX ≡DdK, the first term of Eq.(S21) becomes
α tr[〈(sT dKT g(u) – 〈sT dKT g(u)〉)g(u)sT + (E(u)–E0) Diag[g'(u)] dK ssT〉dKT]NM= α ∑ ∑ [ Dii2 〈g'(si)sj2〉 dKij + Dii Djj dKji Θ[j≤N]
i=1 j=1+ {Dii2 Ωij + 2(Δjj–Δii)〈g'(si)〉Θ[j≤N] – (Dii–Δii)〈g'(si)sj2〉} dKij] dKijNM

= α[ ∑ ∑{ (Dii2–Dii+Δii)〈g'(si)sj2〉 + Dii2 Ωij}dKij2
i=1 j=1
NN+ ∑ ∑{ 2(Δjj–Δii)〈g'(si)〉dKij2 + Dii dKij Djj dKji} ]i=1 j=1

N

M

N

N

≈ ∑ ∑{ aΔii〈g'(si)sj2〉+αΩij}dXij2 + ∑ ∑{ 2(Δjj–Δii)〈g'(si)〉dXij2+αdXijdXji}.
i=1 j=1(S32)

i=1 j=1

In the last line, Dii–2(Dii2–Dii+Δii) ≈ 1 – Dii–1 + Δii ≈ 1 – (1–(a–1)Δii) + Δii = aΔii is used.
In addition, the second and third terms of Eq.(21) become
β{tr[(Φ○(dKT(I,O)))dKT(I,O)] + tr[(Φ○((I,O)TdK))dKT(I,O)]}+ tr[(βΦ○((I,O)T(I,O) – Λ) + γΛ–1)dKTdK]= β{tr[(Φ1○dK1T)dK1T] + tr[(Φ1○dK1,dK2)(dK1,dK2)T]}+ tr[(βΦ1○(I–Λ1) + γΛ1–1)dK1TdK1] + tr[(βΦ2○(–Λ2) + γΛ2–1)dK2TdK2]= β{tr[dK1(Φ1○dK1)] + tr[dK1(Φ1○dK1T) + dK2dK2T]}+ tr[(βΦ1○(I–Λ1) + γΛ1–1)dK1TdK1] + tr[(βΦ2○(–Λ2) + γΛ2–1)dK2TdK2]= βtr[dK1(Φ1○(dK1+dK1T))] + tr[(βΦ1○(I–Λ1) + γΛ1–1)dK1TdK1]+ tr[(β(I–Φ2○Λ2) + γΛ2–1)dK2TdK2]
N

N

N

N

= β ∑ ∑ Φij(dKij+dKji)dKij + ∑ ∑ ( βΦjj(1–Λjj)+γΛjj–1)dKij2
i=1 j=1i=1 j=1

16

NM+ ∑ ∑ ( β(1–ΦjjΛjj)+γΛjj–1)dKij2.

(S33)

i=1 j= N +1

Since Δ1 = (11T+Ω1)†○{βΦ1○(Λ1 – I) – γΛ1–1} = Diag[(1+Ωii)–1 (β Φii (Λii–1) – γΛii–1)]= Diag[(β Φii (Λii–1) – γΛii–1) / a], Eq.(S33) becomes
N

N

N

N

N

M

β ∑ ∑ Φ ij (dKij+dKji)dKij – ∑ ∑ aΔ jj dKij2 + ∑ ∑ ( β(1–ΦjjΛjj)+γΛjj–1)dKij2
i=1 j=1N

i=1 j= N +1

i=1 j=1

NNM≈ ∑ ∑{ (βΦij – aΔjj) dXij2 + βΦij dXji dXij} + ∑ ∑ ( β(1–ΦjjΛjj)+γΛjj–1)dXij2.

(S34)

i=1 j= N +1

i=1 j=1

Therefore, Eq.(S21) becomes
N

M

N

N

d2L ≈ ∑ ∑{ aΔii〈g'(si)sj2〉 + α Ωij}dXij2 + ∑ ∑{ 2(Δjj–Δii)〈g'(si)〉dXij2 + αdXijdXji}
i=1 j=1N

i=1 j=1

NNM+ ∑ ∑{ (βΦij – aΔjj) dXij2 + βΦij dXji dXij} + ∑ ∑ { β(1–ΦjjΛjj)+γΛjj–1}dKij2
i=1 j= N +1

i=1 j=1

N

N

N

N

≈ ∑ ∑{ aΔii〈g'(si)sj2〉+αΩij+2(Δjj–Δii)〈g'(si)〉+βΦij–aΔjj}dXij2 + ∑ ∑ ( α+βΦij)dXijdXji
i=1 j=1i=1 j=1

N

M

N

M

+ ∑ ∑ { aΔii〈g'(si)sj2〉 + αΩij}dXij2 + ∑ ∑ { β(1–ΦjjΛjj)+γΛjj–1}dXij2
i=1 j= N +1N

i=1 j= N +1

N

N

N

= ∑ ∑{ a(Δii〈g'(si)sj2〉–Δjj)+αΩij+2(Δjj–Δii)〈g'(si)〉+βΦij}dXij2 + ∑ ∑ ( α+βΦij)dXijdXji
i=1 j=1N

i=1 j=1M+ ∑ ∑ { (βΦii(Λii–1)–γΛii–1+α
i=1 j= N +1

a−2)〈g'(si)〉 + β(1–ΦjjΛjj)+γΛjj–1}dXij2.
a

(S35)

Therefore, necessary and sufficient conditions to be linearly stable are given by
For 1 ≤ i = j ≤ N,a(Δii(a–1)–Δii) + α(a–1) + α + 2βΦii = aΔii(a–2) + αa + 2βΦii > 0
For 1 ≤ i ≠ j ≤ N,{a(Δii〈g'(si)〉–Δjj) + α〈g'(si)〉 + 2(Δjj–Δii)〈g'(si)〉 + β}
· {a(Δjj〈g'(sj)〉–Δii) + α〈g'(sj)〉 + 2(Δii–Δjj)〈g'(sj)〉 + β} > α+β⟺ α2 〈g'(si)〉 〈g'(sj)〉 + 〈g'(si)〉 {a(Δjj〈g'(sj)〉–Δii) + 2(Δii–Δjj)〈g'(sj)〉 + β}
+ 〈g'(sj)〉 {a(Δii〈g'(si)〉–Δjj) + 2(Δjj–Δii)〈g'(si)〉 + β} > 1⟺ (1–2β)〈g'(si)〉 + a〈g'(si)〉(〈g'(si)〉–1)(Δii+Δjj) + 2β〈g'(si)〉 > 12⟺ (〈g'(si)〉+1)(〈g'(si)〉–1) – 2β〈g'(si)〉(〈g'(si)〉–1) + a〈g'(si)〉(〈g'(si)〉–1)(Δii+Δjj) > 017⟺ (〈g'(si)〉 – 1){〈g'(si)〉 + 1 – 2β〈g'(si)〉 + a〈g'(si)〉(Δii+Δjj)} > 0,
For 1 ≤ i ≤ N, N+1 ≤ j ≤ M,
(βΦii(Λii–1) – γΛii–1 + αa−2)〈g'(si)〉 + β(1 – ΦjjΛjj) + γΛjj–1 > 0.
a

(S36)

The first and second inequalities are satisfied when p0(si) = exp(–b|si|a)/Z. Thus, thelinear stability depends on the third inequality.

S4. Conventional PCA algorithm
Oja’s subspace rule [11] is defined by．
W ∝ 〈u(xT – uTW)〉= K 〈ssT〉 AT – K 〈ssT〉 KTW.

(S37)

When 〈ssT〉 = I, it becomes．
W ∝ K (AT A – KT K) A–1.

(S38)

This is equivalent to the EGHR-β with (α, β, γ) = (0, 1, 0) except positive definite
matrix AAT and kurtosis terms.

18