Chord Function Recognition as Latent State Transition

This paper proposes an unsupervised learning of chord classification aiming at an autonomous recognition of chord functions. In this research, we employ hidden semi-Markov model to incorporate music metrical structure, and in addition, we combine the model with neural network components to embed context information such as beat positions and preceding chord sequences. Experimental results show that the added contexts considerably improve the perplexity. With the help of these neural networks, the proposed model automatically learns hidden states that appropriately represent chord categories. To this purpose, we pre-process the dataset minimally; that is, we only transpose pieces so as not to possess key signatures and ignore octave positions in pitch events. We observe the chord categories effectively cover chords that appeared in the corpus. We further show that the transitions between chord categories reflect the difference of tonalities with a tendency consistent with known chord functions.

However, such textbook theories are too restrictive to represent diverse musical styles. Considering that these theories are derived from the experience of existing music, we are interested in whether a machine can induce inherent regularities in music by statistical learning. Unlike direct applications as music generation systems [11,12], dataoriented, statistical acquisition of chord functions has been rarely investigated [13][14][15][16][17]. It has been reported that the statistically induced regularities considerably agree with known textbook theories, and it is expected to find finegrained patterns that are characteristic of particular genres or composers [13][14][15][16][17].
These previous works, however, required the chord labeling or grouping as a pre-process, which includes chord segmentation [13,16], modulation segmentation [14,16], and chord annotation with scale degrees or Berklee chords [14,15,17]. Although harmony analyses in textbooks are based on given chord names, the correspondence between chord names and surface pitch events in scores is not trivial and diverse depending on musical styles. Thus, obtaining consistent Berklee annotations [18] or key assignments [19] were found to be more difficult than generally expected. In addition, pre-processed modulation detection or chord degree assignment ignores the mutual dependency of chord functions and keys. For example, G major chord is dominant in the C major key, while it becomes tonic in the G major key. Thus, the same chord may possess different regularities of progressions according to tonality [2]. Based on this knowledge, previous works for chord function identification have simplified the problem by removing the dependency for keys in advance. However, we regard that these functional changes of codes are also important features of music that should be captured in data-oriented statistical models. This article is part of the topical collection "Evolutionary Art and Music" guest edited by Aniko Ekart, Juan Romero and Tiago Martins.
We present a model that automatically learns chord categories and segments surface notes to be classified into the categories without relying on pre-defined chord symbols and key-assignments by utilizing an extension of hidden semi-Markov models (HSMM), where the chord categories are learned as hidden states that are used to predict the next coming chord in a chord sequence. As textbook theory required the concept of chord functions to explain chord progressions, we consider the characteristics of chords and the progressions between chords to be important musical features worth learning from the data.
Although our future goal is to apply the data-driven model to a variety of compositions to reveal the characteristics of individual composers, we start with J. S. Bach's four-part chorales in this study. In this research, we avoid any heuristic pre-processing and only apply ones that are simply based on musical notation; in particular, we transpose pieces so as not to possess a key signature 1 and ignore octave positions.
Technically, the task is a kind of sequential tagging, and thus it is close to the unsupervised part-of-speech (POS) tag induction in the natural language processing (NLP). We employ the unsupervised neural HMM, which is known to be one of the prominent models proposed for the POS tag induction [20]. The most advantageous feature of the neural HMM is that we can easily embed additional contexts. We customize the model and its additional contexts to be suitable for the music. Furthermore, to make the model be in accordance with metrical structures, we extend it to the neural hidden semi-Markov model (HSMM) 2 .
We use perplexity as the evaluation metric based on previous works [15,17]. Experiments show that our model appropriately segments and classifies surface pitch-classes, especially with the smallest perplexity. Additional contexts with neural network modeling are shown to be effective to improve the perplexity. In addition, we show the transitions between categories reflect the difference of tonalities when we count them by separating pieces into groups of majors, minors, and dorian scales. This paper is organized as follows. In Sect. "Related Work", we review related studies. We introduce the proposed model in Sect. "Unsupervised Neural Hidden Semi-Markov Model for ChordClassification". Then, we show the experimental results in Sect. "Experiments". Finally, we summarize our contributions in Sect. "Conclusion".

Related Work
Unsupervised learning of harmony or chord progressions has been studied through clustering methods [13,14], HMMs [15][16][17], and Probabilistic Context-free Grammars (PCFG) [15]. These studies aimed at acquiring the progression rules, not based on possibly subjective human annotation, but relying on data-driven analysis.
They have found similarities between statistically induced clusters and chord functions in textbook theories, when compared models with the same number of clusters as the textbook [14][15][16]. Especially in HMM based models, not only obtained clusters but also their state transition property was found to have similarity to the known chord functions [16,17]. However, selecting an optimal number of clusters would not be trivial. White and Quinn [16] proposed a methodology to find the best number of states for HMMs. They adopted the k-medoids clustering over hundreds of HMMs with different initial parameterizations. And then, they assumed that the number of hidden state N that showed higher silhouette width (i.e., clearer cluster boundary) than N − 1 and N + 1 was appropriate. Uehara et al. [17] proposed another tandem approach to find the number of hidden states and detect modulation segments. On the other hand, Jacoby et al. [14] investigated the balance of the model by introducing the optimal complexity-accuracy curve, rather than fixing a particular number of hidden states. Similarly, Tsushima et al. [15] evaluated the generative models by the perplexity, which is a common metric for estimating generalization performance of a statistical model, and found that the larger number of hidden states leads a better score. Based on the latter approaches, we train our model with multiple numbers of hidden states and evaluate them by perplexity.
Preceding studies mentioned focused on harmonic structure solely. In other words, the duration of chords was ignored. However, a recent work proposed a supervised learning of the combined model of harmony and rhythm, and reported the efficacy of rhythmical information [21]. We also incorporate metrical information into our unsupervised learning by a semi-Markov model. Not only employing a semi-Markov model, we adopt an extension with neural networks based on the neural HMM. An important strength of the neural HMM is the seamless integration with additional contexts. Tran et al. [20] introduced the two additional contexts: embedded feature of preceding observations by the Long-Short Term Memory (LSTM) [27] and morphological information via character Convolutional Neural Networks. Despite its simple framework, the neural HMM outperformed even the highly polished Bayesian mixture model [22], or the hierarchical Pitman-Yor Process HMM [23]. Additional contexts in the neural HMM are also can be used for the unsupervised neural HSMM in the same manner.

SN Computer Science
Unlike the Hybrid DNN(Deep Neural Network)-HMM [24] that converts a pre-trained GMM(Gaussian Mixture Model)-HMM to a DNN-HMM or the Tandem model that combines features by a supervised DNN [25] to an HMM, the neural HMM is seamless and fully unsupervised. Our extension, i.e., the unsupervised neural HSMM differs from the Recurrent HSMM that makes use of the bi-LSTM for the purpose of reducing the error in the variational approximation [26].

Hidden Semi-Markov Model
We aim at presenting a model that predicts chord segments, chord categories and chord progressions simultaneously. To achieve this, we employ hidden semi-Markov model (HSMM), which is an extension of hidden Markov model (HMM), where lies a Markov chain of hidden states behind an observable sequence. The notion of the duration of a hidden state is introduced by the "semi-Markov" extension [28,[31][32][33]. While there can be various ways to model the duration of a state, usually one of the following three models is used for computational efficiency [34]: Explicit duration HMM [28,31], Variable transition HMM [32], and Residential-time HMM [33]. We select the Residentialtime HMM [33,34] 3 since the computational complexity of which is the smallest of the three 4 .
Residential-time HMM assumes that a hidden state transition is independent of the duration of the previous hidden state and the duration of a hidden state is a discrete random variable. With this assumption, hidden state transition is described as follows 5 where i, j ∈ {0, 1, ..., S} are state indices, S is the number of hidden states, ∈ {1, 2, … , D} is the discrete duration time, and D is the maximum duration of a hidden state. We call a particular value of i or j as a hidden state index. While a hidden state index itself is just a serial number in hidden states of H(S)MM, a particular hidden state is expected to be a chord category after the training. Furthermore, a i,(j, � ) can be decomposed into transition probability and duration probability as follows.
Note that a hidden state changes to another only when the remaining duration = 1 . Therefore, the duration probability determines a hidden state duration, and the self-transition probability a ii is always zero.
As is shown in the graphical representation of HSMM (Fig. 1a), each hidden state ( z t] = i ) can produce multiple observations. Then, when a hidden state changes to another one (j), the model calculates the duration probability of the next state p j as well as the transition probability a ij . However, there could be multiple possible combinations of the hidden state and the duration. In Fig. 1b, the number of hidden states is 3 and the maximum duration of a hidden state is 2; when the hidden state index is i = 0 and (remaining) duration is = 1 at the time step t, possible conditions for the previous time step t − 1 are: (i = 0, = 2) , (i = 1, = 1) , and (i = 2, = 1) . Note that if the current remaining duration is = 1 , the hidden state must change to another one (dashed lines) at the next time step. On the other hand, if > 1 , the hidden state continues into the next time step (solid lines).

Framework
We expect that hidden states represent chord categories, which are not given a priori and thus learned in an unsupervised manner. With the Markov property assumption, hidden states (chord categories) are learned to govern the progression of the forthcoming chord category. We argue that an appropriate set of chord symbols would be diverse along with targeted pieces and thus adopt unsupervised learning.
We give a more detailed description of the proposed framework by showing an example in Fig. 2. Since HMM (HSMM) takes discrete time series, we set a time step by every sixteenth note, which is the minimum duration of note in the corpus, except for quite few exceptions of 32nd notes. We obtain a pitch-class vector for each segment 6 . For example, when the time step is at 25, the pitch names of notes contained in the segment are (C, D, G), and thus corresponding pitch class vector is (0, 2, 7).
We build a table of vocabulary of observations, and call them i.e., tokens, from every combination of the four-part pitch-classes. For example, the token index for (0, 2, 7) is set to 16. The whole vocabulary is shown in Appendix a i,(j, � ) = a ij p j � 3 While Residential-time HMM was originally proposed as an efficient approximation of Explicit duration HMM [33], the same author named it the "Residential-time HMM" in the latter article [34]. 4 The computational complexities of the three models are follows [33]. ( Table 4). As are shown in the example, observations may include passing tones. We expect that a hidden state works as a chord category behind these raw pitch-class vectors. Looking at time step 25 again, we can see that the index of the hidden state is 7. Here, the hidden state is lasting from time step 25 to 28, which includes three types of pitch class vectors, (0, 2, 7), (2,7,11) and (2,5,7,11), but all the three can be interpreted as a part of G major chord.
A hidden state emits a token with probability b ik , where i is a hidden state index and k is a token index. In the case of time step 25, pitch-class vector (0, 2, 7) (the token index of which is 16) is emitted with the probability of b 7,16 , from the hidden state indexed 7. A hidden state is changed to another one with the probability of transition probability a ij multiplied by duration probability p j , where i is the index of preceding hidden state, j and are the index and duration of the forthcoming hidden state respectively. The probability of the transition from hidden state 7 to hidden state 3 at time step 28 to 29 is then a 7,3 p 3,8 .
Note again that hidden states and their duration are not given a priori, but obtained by unsupervised learning. More precisely, we first optimize parameters that determine transition, duration, and emission categorical distributions so that a tuned model gives higher likelihoods to targeted pieces, and then obtain the best hidden state sequence based on the model. Next, we give a detailed description of the architecture of neural networks used in the neural HSMM in the following section.

Architecture of Neural Hidden Semi-Markov Model
Different from the conventional HSMM that has categorical parameters just as matrices of parameters, neural HSMM equips neural network components as functions for calculating transition, duration, and emission categorical distributions. As mentioned in the previous section, the role of each distribution is as follows.

SN Computer Science
-Transition distribution: the probability of transition from a hidden state to another one, which corresponds to transition of chord categories. -Duration distribution: the probability of duration of a hidden state, which corresponds to duration of a chord category. -Emission distribution: the probability for a hidden state to emit a particular token. In our case, this is the probability for a chord category to emit a particular pitchclass vector.
Besides these three distributions, there is a special case of transition distribution i.e., initial state distribution, described in Sect. "Initial Hidden State Probability". The same graphical representation ( Fig. 1a) applies to neural HSMMs, but the categorical distributions are obtained as outputs of neural networks that can employ additional musical contexts for the calculation. We show these networks and additional contexts in the following paragraphs: hidden state transition probability ("Hidden State Transition Probability), initial state probability ("Initial Hidden State Probability"), duration probability ("Duration Probability"), and emission probability ("Emission Probability") 7 .

Hidden State Transition Probability
We denote a hidden state at time step t as z t and state index of which as i or j. Given the number of hidden state S, each hidden state i ∈ S has a transition distribution to the next hidden state j ∈ S , which is denoted as a ij , and thus ∑ j a ij = 1 . Note that, in HSMM (Residential-time HMM), the self-transition a ii is set to zero, since not transition distribution but duration distribution manages the duration of a hidden state. Instead of having a ij just as a single learnable value, we prepare a function that calculates it by employing neural network components as follows.
where MLP 3 is a 3-layer Multi Layer Perceptron that has one input layer, two hidden layers and one output layer with a hyperbolic tangent (tanh ) activation function after each hidden layer 8 . The output layer of this function is softmax 9 , to satisfy the condition ∑ j a ij = 1 . The output layer size of the transition MLP 2 in (1) should be one smaller than the number of states ( S − 1 ), since it does not allow self-transitions.
(1)   Table 4 in Appendix. Residential time represents a remaining duration of a hidden state i, which decreases by 1 in the same hid-den state. Time steps increment by the sixteenth-note width segment. The notation z [13∶16] = 2 represents that the hidden state index is i = 2 on time steps 13-16. a ij is a transition probability and a 2,3 represent that of hidden state index 2 to 3. b ik is an emission probability and b 2,12 is that of token index 12 emitted from hidden state index 2. p j is a duration probability and p 3,2 represents that the next hidden state index is 3 and the duration of which is 2 7 Summary of notations (Table 5) and model settings (Table 6) can be found in Appendix. 8 Similarly, MLP 2 denotes a MLP with one hidden layer, and we use the same notation in following sections as well.
Since a ij corresponds to the transition probability of hidden state i to j, the input for the network includes a feature for hidden state i. We provide a hidden state embedding s i , which is a learnable vector associated to hidden state index i. This hidden state embedding is also jointly used in the network for duration and emission probability (described in Sects. "Duration Probability" and "Emission Probability"), and would help to capture relationships between hidden states (chord categories) and emissions (observed pitchclass vectors). Taking further advantage of neural HMM, we introduce two types of additional contexts for the calculation of hidden state transition probability, i.e., a pitch-class histogram and an embedded feature of preceding pitch-class vectors, and feed a concatenation of the two contexts c t to (1).
Additional context 1: Pitch-class histogram ( ) The first additional context is a pitch-class histogram that is the frequency of occurrence of each pitch-classes calculated from an entire phrase, which represents tonality. We split a piece into chord sequences at each fermata (point d'orgue) that works as a full-stop marker of a lyric, and then sum up the duration of each pitch-class to obtain a pitch-class histogram of a sequence. Even though the main tonality of a piece would be easily distinguished by the key signature 10 , there would be local modulations and thus we employ such a histogram as an indirect information of local tonality. We obtain the feature of a pitchclass histogram of a chord sequence r histo as follows.
where v histo ∈ ℝ 12 is the raw pitch-class histogram of a sequence.
(2) r histo = MLP 2 (v histo ) Additional context 2: Embedded feature of preceding observations by the LSTM ( ) The second additional context is an embedded feature of preceding observations by the Long-Short Term Memory (LSTM) [27], based on the idea from [20]. LSTM is a variant of Recurrent Neural Network (RNN) that recursively feeds observations, and thus it is suitable for a feature function for sequential data. In our case, it corresponds to a feature of a temporal sequence (from t = 1 to current time step t) of pitch-class vectors. At each time step, we input an observation embedding o k that associated with the observation v pitch k at t 11 to LSTM as follows (3)(4).
In calculation of an observation embedding (3), v pitch k is a binary pitch-class vector, that is a 12-dimensional vector consisting of 1/0. For example, if the pitch-class vector is (2,5,7,11), the corresponding binary pitch-class vector is (0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1). We expect that binary pitch-class vectors give raw information of observations, instead of symbolized token indices, and it is also used in emission distribution network described in Sect. "Emission Probability". Finally, we concatenate the two contexts and obtain the context vector as (5), where [; ] denotes vector concatenation. This context vector is applied to the transition probability network (1).
We illustrate the network for calculating transition probability a ij in Fig. 3, which we have described in this section. Note that transition probability thus dynamically varies by  11 In this example, x t = k , which means that the token index for x t is k. 10 However, the main tonalities of pieces that are written in church modes (e.g., dorian) would not correspond to the key signature system of modern tonality. In addition, major/minor is not distinguished by the key signature solely. Furthermore, some minor pieces retain the feature of church modes [3] and frequently proceed relative major keys [13,17].

SN Computer Science
contexts in neural HSMM, though it is static in conventional models. Therefore, while we train a model without distinguishing major/minor pieces or local modulations, we expect that the additional contexts help to automatically adjust transition probabilities.

Initial Hidden State Probability
Since the first pitch-class vector at time step 1 does not have a preceding observation, this special distribution for the initial hidden state is provided. We set an external initial hidden state z 0 at time step 0 that has no emission and duration 1, to meet the initial boundary condition of HSMM. The initial hidden state probability i is calculated as follows. This is the only component remained almost the same as the conventional model except that it is normalized by the sof tmax function.
where i is a hidden state index, ∈ ℝ S is the learnable weight vector, and S is the number of hidden states.

Duration Probability
The duration probability determines how long the hidden state i (i.e., a chord category) will reside in the same hidden state. The network for calculating it is as follows. The output layer size of the duration MLP 2 in (6) is 16 that is a maximum duration length that corresponds to a whole note. The hidden state embedding s i that is used in the network for transition probability is used here again, and is expected to help to jointly learn the transition and duration distribution of hidden state i.

Emission Probability
We employ a discrete HSMM in this study, therefore, each observation is associated with token index k ∈ V in the vocabulary as mentioned in the example in Sect. "Framework". The probability that hidden state i yields an token index k (b ik ) is calculated as follows, based on [20].
Additional context 4: Observation embedding from binary pitch-class vector ( ) We use the same observation embedding o k as one used in the context for transition probability (3) for the calculation of emission probability (8). Instead of employing weight matrices or MLP s, the emission probability is learned as the dot product of a hidden state embedding and an observation embedding. The bias value l k ∈ ℝ would help to consider the frequency of each chord. In contrast to the conventional discrete HMM that cannot use the raw information about the observation once converted into a token index, our model does not lose it by the observation embedding.

Optimization
It is known that there is no analytical solution to obtain an optimal parameterization of an HMM that maximizes the likelihood of the observed sequence [28]. Although the widely used Baum-Welch re-estimation algorithm (that is a kind of the expectation-maximization algorithm) is able to obtain a locally optimal parameterization, it would be stuck in bad local optima when the likelihood surface is complex [28]. On the other hand, we can utilize gradient based methods to maximize the likelihood or the marginal probability 1∶T is an observed sequence since it is an optimization problem [28][29][30]. We employ the known dynamic programming named the forward algorithm for H(S)MMs to calculate the marginal probability [28,34].
Neural networks models are generally trained by gradient based optimizers as well as the backpropagation for the gradient computing. Recently, along with the success of neural networks, efficient optimizers also have been proposed. By implementing the HSMM as a neural network, we can naturally utilize the gradient based optimization [29]. In this study, we use one of the latest optimizer RAdam [35], which adjusts the learning rate automatically.

Forward Algorithm for Hidden Semi-Markov Models
To conduct the gradient based optimization, we calculate the marginal probability of an observed sequence ln P( 1∶T ) by the forward algorithm, which marginalizes the possible paths (as previously illustrated in Fig. 1b) of states and durations. We give the forward algorithm for the Residentialtime HMM in the following.
To explain the duration of hidden states, we borrow the notation used by Yu et al. [34].
:a hidden state lasts at latest from t 1 and ends at time t 2 .
Among variants of HSMMs, the Residential-Time HMM does not allow the self-transition, and a transition to another hidden state is required when the "remaining" duration is = 1 . On the other hand, when > 1 , the duration just decrements at each time step.
The t (j, ) represents the forward probability that the hidden state at t is j and the remaining duration of which is , and output tokens from t = 1 to t, i.e., 1∶t are observed. It is decomposed as follows.
where P( |z t = j) = p j is duration probability 12 , P(x t = k|z t = j) = b jk is emission probability, and P(z t = j|z t−1 = i) = a ij is transition probability. Note that when the state at t is j, there are two possibilities: (i) the hidden state at time t − 1 is also j with the remaining duration at there is + 1 , (ii) the hidden state at time t − 1 is i (i ≠ j) and the remaining duration is 1 then is transferred another hidden state j at the time t.
We apply scaling by replacing t (i, ) to the conditional probability, similar to HMM [34], and obtain a modified forward algorithm.
As can be seen from (9), the marginal probability is obtained by P where T is the length of an entire observation sequence, which means we can train the model only by executing the forward algorithm, without the backward algorithm usually used in EM algorithm. Note that C t is obtained by summing up (10) about all j and , since ∑ j ∑̂t (j, ) = 1 . Unlike HMMs, we set the initial probability to yield the initial hidden state z 0 that does not yield observation by considering the initial boundary condition: 0 (i, 1) = P(z 0 = i) and 0 (i, ) = 0 for > 1 [34]. Therefore, 0 =̂0 and C 0 = 1.0 at t = 0.

Experiments
In this section, we show the experimental results. We first describe the dataset (Sect. "Dataset") and experimental setups (Sect. "Experimental Setups"). Thereafter, we show the evaluation result by perplexity in Sect. "Evaluation by Perplexity". We give qualitative analysis for obtained chord categories and their progressions in Sect. "Qualitative Analysis for Induced Clusters". Finally, we show a detailed exemplification and discussion of BWV267 with reference of a human annotation in Sect. "Discussion on an Analysis by the Model".

Dataset
We use J.S.Bach's four-part chorales by the Riemenschneider numbering system (1-371) from the Music21 Corpus [37] as our dataset resource. At first, we have removed duplicated pieces according to the analysis by Dahn [38]. We only ∈ D is a discrete value of remaining time steps, and D is the maximum duration.

SN Computer Science
contain four-four time (4/4) pieces 13 , thus the number of pieces is 290 (Table 1). We regard each phrase as an independent sequence, however, when we provide these phrases to training, evaluation, and testing set, we give randomness over pieces (not over phrases), since one piece may contain similar phrases.

Experimental Setups
In this research, we disregard the difference of absolute pitch so that we transpose all the major pieces to C major and all the minor ones a minor. Some pieces of the corpus are written in dorian scale without key signature. From the viewpoint of modern tonality, we regard them in d minor; however, we do not shift them to a minor. Also, local modulations are shifted in a similar way, so as to preserve the relative positions of constituent notes. Finally, we ignore the difference by octave in pitch events and inversion in chords.
We assign token indices (k) for pitch-classes the accumulated pitch-duration of them in the dataset over 95%, and then, remained chords are merged to " ". The obtained vocabulary size was 80, including " ". In training, we set the mini-batch 8. We train models up to 500 epochs, however, the process is stopped when the lowest loss on the evaluation set is not updated over 20 epochs. We apply the dropout [36] to each hidden layer of MLPs, which randomly ignores some ratio of neurons (12.5% in our case) during training to avoid over-fitting.

Evaluation Metric
Following previous works [15,17], we use perplexity as an evaluation metric, defined by the following equation.
where (x 1 , … x T ) is a sequence of output tokens and is a set of model parameters. The smaller perplexity for test data means the better generalization performance since perplexity corresponds to log-average of inverse probability, and thus it is commonly used for evaluating probabilistic models. We examine the performance on multiple numbers of hidden states: 3 to 16. We conduct the experiments with three different random seeds of {0, 1, 2} for each number of hidden states, and report the averaged scores in Sect. "Evaluation by Perplexity".

Perplexity Scores
In addition to the neural HSMM, we implement a baseline model that represents probabilities as simply learnable weight vectors or matrices with softmax output layers. The baseline model is almost "Non-Neural" but an HSMM tuned by the same gradient-based optimizer. We compare the proposed model with this baseline to see how efficiently the elaborated neural network components work. In Fig. 4, we can see that the elaborated neural models considerably outperformed the baselines. The ablation study also shows the efficacy of additional contexts as shown in Table 2  the perplexity. Among them, removing 14 the observation representation by pitch-classes ( − ) led a significant drop. Since we do not employ MLPs for the architecture of calculation of emission probability but force a model to learn relationships between hidden states and vocabularies directly, the information of raw pitch-classes would help the learning.  (1)

Qualitative Analysis for Induced Clusters
In this section, we give qualitative analysis and discussion about the induced chord clusters. We focus on the eight-state models since we have transposed pieces not to have key signatures and have assumed that each cluster would roughly correspond to a triad on the diatonic scale, or otherwise to the rests.

Induced Clusters and Model's Perplexities
We investigate if we can find clearer clusters in the higher scored model in obtained examples. We show the emission probabilities of the best neural HSMM, the worst neural HSMM, and the best baseline model in Fig. 5. Note that we have executed three trials of training for each model with a different random seed among {0, 1, 2}. Here the "best" model means that the evaluation perplexity of which is the smallest among the three trials. We observe that the clusters obtained by the best scored neural HMM mainly consist of the chords on diatonic scales (C major and a minor) and its commonlyused borrowed chords such as D major chord. We show the emission probabilities of the top three emissions for each hidden state index in Fig. 5 and summarize them in Table 3. For ease of readability, we name chord categories by the name of the most frequent output chord for each hidden state index. Although the most frequent chord is representative of the category, it does not correspond to a unique chord but a set of emission probabilities; therefore, we add "hat" to it. The cluster Â may emerge from two reasons: the common use of the Picardy third, or the dominant of dorian mode pieces. It is worth noting that seventh chords and passing chords are merged into appropriate clusters, e.g., C.D.E.G in the same hidden state for C.E.G (state3), and D.F.G.B for D.G.B (state7).
Even though employing the same neural HSMM, the worst perplexity model (in the middle of Fig. 5) seems to be less appropriate than the best model. For example, C.E.A is mixed up with C.F.A and E.G♯(A♭).B in state7, and C.E.G (state4) and passing chords around it (state0) are separated. Although we would be able to choose a model by perplexity, we admit that difficulty for obtaining the global optimum still exists even though we employ the efficient gradient based optimizer. The best baseline model (the evaluation perplexity of which is 9.92) is still worse than the worst neural HSMM (9.18). Not only the best baseline model scored a worse perplexity, but it possessed more miscellaneous clusters than neural HSMMs. However, we can observe that C.E.G and D.F.G.B of state3, or D.G.B and D.F♯(G♭).A of state5, in the same tonality, would tend to be merged into the same category in the baseline HSMM.  Table 3 Hidden State Transitions Since chord functions lie in the regularity of chord progressions, we try to find them by investigating the transitions between obtained chord categories. We count the hidden state transitions on major pieces, minor pieces and dorian pieces separately to see whether the model appropriately reflects the difference of tonalities in its state transitions. Note that, instead of examining the hidden state transition probability directly, we count the number of transitions after decoding the sequence of hidden states by the Viterbi algorithm since the hidden state transition probability changes by contexts.
We show the result of hidden state transition properties by the best scored neural HSMM in Fig. 6 (the emission probability of which is shown at the top of Fig. 5). In major pieces, the tendency that {state4:F, state2:d} (subdominant) → state7:Ĝ (dominant) → state3:Ĉ (tonic) is noticeable. While in minor pieces, the tendency described above suggests existence of the relative major keys, which is consistent with previous works [13,17]. In addition, a strong transition from state5:Ê (dominant) → state6:â (tonic) are observed. Unlike in major pieces, the transition from state2:d̂ (subdominant)

Discussion on an Analysis by the Model
In this section, we discuss the adequacy of our model from an obtained result. We select BWV267 from our testing set, since a human analysis on which is publicly available [37] 15 . Although our aim is not to reproduce the human analysis, we consult it as a possible interpretation for local modulations and chord annotations. Note that we normalized the score to have no key signature, and thus the main key became C major, while the original key was G major. We show extractions of the result in Fig. 7 16 .
According to the human analysis [37], there are local modulations to {F major, G major, d minor, g minor} keys in this piece.
In this example, we could see that the model found an effective set of clusters to cover chords appeared in a piece including local modulations. For example, in the section of G major key (phrase No.4, time step 21-32), we observed the cluster of D̂ (state0), which was a borrowed chord seen from C major key. This cluster D̂ was used again in the section of g minor key (phrase No.7, time step 7-32). Interestingly, g minor chords, i.e., D.G.A♯ (B♭ ) ( k = 21 ) and D.F.G.A♯(B♭ ) ( k = 60 ), were classified into the same cluster as Ĝ (state7), since g minor chord and G major one shared the dominant (state0: D). Similarly, in the section of F major (phrase No.4, time step 1-20), C major dominant seventh chords, i.e., C.E.G.A♯ (B♭ ) ( k = 24 ), were classified into the same cluster as C major chords (state3: Ĉ). According to Schoenberg, local modulations basically remain in closely related keys and thus can be analyzed by altered chords [3]. In the case of g/G chords and C/C7 chords described above, these alternations (g or C7) do not change the probable progression of chords, and thus would be classified into the same function as the one on the diatonic scale (G or C).
In the longest modulation section of d minor (phrase No.6), we could see that the cluster of tonic chord (state3: Ĉ) has not appeared and that the cluster of d minor chord (state2:d) tended to proceed to A major chords (state1:Â).
Since we ignored inversions of chords in this study, we basically could not capture the difference between them. We admit that it needs further consideration, and left it as a future work. In addition, we observed more " " tokens ( k = 79 ) in modulation sections that hindered the analysis. " " or " " token has been introduced in discrete HMMs to avoid useless increase of dimension of the emission distribution. However, once an observation is converted into " ", the raw information about it (which is pitch-class vector in our case) will be lost. More effective treatment for less frequent observations or more direct representation for emission distribution is another important future work for us.

Conclusion
In this paper, we aimed to obtain the regularity of chord progressions by the data-oriented unsupervised learning. Chord functions were introduced to explain typical progressions a posteriori, and thus music is not constrained by such functions. Thus we have been interested in finding statistically plausible chord functions by machines.
Harmony analysis consists of the following four processes; (i) defining an appropriate set of chord labels (e.g., chord degrees or Berklee chords in traditional approaches), (ii) determining the scopes of chords and their labels from surface structures of a piece, (iii) identifying the sequence of chords based on the determined labels, (iv) and finally, analyzing the labeled chord sequence by examining transitions between chords, i.e., chord functions. Our objective is to let computers execute these processes, without relying on human annotations or pre-defined labels. In our experiments, we have only pre-processed the source scores by shifting keys either to C major or to a minor, and ignored the difference by octave in pitch events. The pre-process could be executed automatically without relying on external music feature extraction tools.
We regarded that the states in hidden Markov model (HMM) represented chord categories, and we employed an extended HMM with the notion of duration of a hidden state, called hidden semi-Markov model (HSMM). Furthermore, we utilized neural networks for calculating categorical distributions to embed the pitch-class distribution, preceding chord sequences, beat information, and so on, as additional contexts.
The experimental results satisfactorily showed the efficacy of such contexts in terms of perplexity in comparison with the baseline and other ablation studies (Sect. "Evaluation by Perplexity"). We successfully obtained clear chord clusters with the best scored proposed model where seventh chords and borrowed chords were classified into appropriate clusters (Sect. "Induced Clusters and Model's Perplexities"). Also, we have shown a full exemplification of BWV267 (Sect. "Discussion on an Analysis by the Model").
As was mentioned by Hugo Riemann [2], "In the change of these function (tonic, dominant, and sub-dominant) lies the essence of modulation". According to this, the distinction of keys should not be given prior to the chord analysis, but the modulations should be externalized as a result of the analysis. We followed this aphorism and have avoided the independent detection of modulation. Instead, we have characterized each hidden state, counting the number of their transitions, and dividing the experimental results into majors, minors, and dorian scale. In this regard, we have not automatically detected the keys, however, we have found that the distribution of state transitions were consistent with the known classical theory, in an unsupervised manner (Sect. "Hidden State Transitions").
We plan to improve the proposed methodology to incorporate modulation detection, to integrate with the chord analysis, as our future work. Also, we intend to apply our model to pieces that may not be faithful to classical musicology. We need further extensions to overcome the out-of-vocabulary issue described in Sect. "Discussion on an Analysis by the Model" for pieces that accompany more elaborations than four-part chorales. In addition, there are usually few pieces in a particular genre by a specific composer. In such a case, learned parameters of transition probabilities and emission probabilities are expected to be good initial values (pre-trained information) for fine-tuning, which is known to be effective when a sufficient amount of training data is not available. For this reason, we used Bach's four-part chorales that contain a substantial amount of coherently structured pieces composed at the beginning of the tonal music. We hope the future work of extensive music analysis by machines would contribute to expanding the scope of discovering cultural evolution in music.

A Appendices
See Tables 4, 5 and 6.   Table 6 The size of layers and related equations Layer Size Eq.