Compact representation of temporal processes in echosounder time series via matrix decomposition

Echosounders are high-frequency sonar systems widely used to observe mid-trophic level animals in the ocean. The recent deluge of echosounder data from diverse ocean observing platforms has created unprecedented opportunities to study the marine ecosystems at broad scales. However, there is a critical lack of methods capable of automatic and adaptive extraction of ecologically relevant spatio-temporal structures from echosounder observation, limiting effective and wider use of these rich datasets in marine ecological research. Here we present a data-driven methodology based on matrix decomposition that builds a compact representation of long-term echosounder time series using intrinsic features in the data, and demonstrate its utility by analyzing an example multi-frequency dataset from the northeast Pacific Ocean. We show that Principal Component Pursuit (PCP) successfully removes noise interference from the data, and that a temporally smooth Nonnegative Matrix Factorization (tsNMF) automatically discovers a small number of distinct daily echogram patterns, whose time-varying linear combination (activation) reconstructs the dominant structures in the original time series. This low-rank representation is more tractable and interpretable than the original time series. It is also suitable for visualization and systematic analysis with other ocean variables such as currents. Unlike existing echo analysis methods that rely on fixed, handcrafted rules, the data-driven and thus adaptable nature of our methodology is well-suited for analyzing data collected from unfamiliar ecosystems or ecosystems undergoing rapid changes in the changing climate. Future developments and applications based on this work will catalyze advancements in marine ecology by providing robust time series analytics for large-scale, acoustics-based biological observation in the ocean.

diverse ocean observing platforms has created unprecedented opportunities to study the marine ecosystems at broad scales. However, there is a critical lack of methods capable of automatic and adaptive extraction of ecologically relevant spatio-temporal structures from echosounder observation, thwarting effective use of these rich data in marine ecological research.
2. Here we present a data-driven methodology based on matrix decomposition that builds a compact representation of long-term echosounder time series using intrinsic features in the data. We use Principal Component Pursuit (PCP) to remove noisy outliers and employ temporally smooth Nonnegative Matrix Factorization (tsNMF) to discover high-level spatio-temporal structures, such as scattering layers and diel migration patterns, from the data. We demonstrate the utility of this methodology by analyzing a multi-frequency time series from an echosounder moored at the shelf break of the northeast Pacific Ocean.
3. Our methodology successfully removes noise interference and automatically discovers a small number of distinct daily echogram patterns, whose time-varying linear combination (activation) reconstructs the dominant structures in the original time series. Each pattern represents a distribution of echo energy across time (hour of day) and space (depth) that co-varies throughout the observa-2 Compact representation of echosounder time series tion period, with multi-frequency features suggestive of scatterer identity. The pattern activations further offer quantitative summarization of temporal processes embedded in the data, in a format that is suitable for visualization and systematic analysis with other ocean variables such as currents.
4. We show that matrix decomposition methods are powerful techniques that could transform complex echo observation into low-dimensional components that are more tractable and interpretable than the original data. Unlike existing echo analysis methods that rely on fixed, handcrafted rules, the data-driven and thus adaptable nature of our methodology is well-suited for analyzing data collected from unfamiliar ecosystems or ecosystems undergoing rapid changes in the changing climate. Future developments and applications based on this work will catalyze advancements in marine ecology by providing robust time series analytics for large-scale, acoustics-based biological observation in the ocean.

I. INTRODUCTION
Sound is extensively used to study life in the ocean (Medwin and Clay, 1998). Compared to net-based sampling, which can only be conducted at discrete times and locations, echosounders boast the capability to "connect the dots" by delivering continuous active acoustic observation across time and space. This advantage has made them an indispensable tool in modern ecological and fisheries studies, in particular for collecting information about mid-trophic level organisms that are otherwise difficult to observe effectively at large scale (Benoit-Bird and Lawson, 2016;Handegard et al., 2013).
Technological advancements have resulted in a deluge of echosounder data in the past decade, thanks to the rapid and successful integration of autonomous echosounders with a wide variety of ocean observing platforms, including moorings and autonomous surface and underwater vehicles (De Robertis et al., 2018;Dunlop et al., 2018;Greene et al., 2014;Moline et al., 2015;Suberg et al., 2014). The significantly broader spatial and temporal coverage of the datasets provides an unprecedented opportunity to study the response of the marine ecosystems to climate variability at scales never before possible. An example is the network of six moored echosounders maintained by the Ocean Observatories Initiative (OOI) along the coast of Northeast Pacific, which have been delivering data continuously since 2015 (Fig. 1A).
Another example is the experimental fisheries survey carried out by a fleet of autonomous surface vehicles sampling the west coast of the U.S. and Canada simultaneously in the summer of 2018-2019 (Chu et al., 2019). However, the massive volume and complexity of 4 Compact representation of echosounder time series these new data have overwhelmed the conventional echosounder data analysis pipelines that rely heavily on manual processing, thwarting rapid progress in marine ecological research.
A critical need in the analysis of large echosounder datasets is the automatic extraction of biological information, such as organism identity, biomass, and high-level spatio-temporal distribution patterns of organisms, without excessive human intervention. Currently, the majority of echo data processing procedures require visual inspection and labeling of echograms (image formed by echo returns, e.g., Fig. 2A) by human experts (Fleischer et al., 2012), prompting not only scalability issues with increasing data volume but also subjective bias and reproducibility challenges.
Many echo classification methods were proposed to automate the procedure of discerning the taxonomic identity of acoustically observed biological aggregations. A key feature commonly exploited is the relative strength of echo across sonar frequency, which varies significantly as a function of the size, shape, orientation, and material properties of marine animals (Benoit-Bird and Lawson, 2016;Medwin and Clay, 1998). These include classification rules derived using empirical echo measurements and physics-based scattering models, such as the popular frequency-differencing and related Bayesian inference methods (e.g., De Robertis et al., 2010;Jech and Michaels, 2006;Korneliussen et al., 2016), as well as supervised and semi-supervised machine learning methods, such as random forests and neural networks (e.g., Brautaset et al., 2020;Fallon et al., 2016;Haralabous and Georgakarakos, 1996;Woillez et al., 2012). However, these classification methods can be ill-fitted for largescale echosounder data collected on autonomous platforms, due to challenges associated with collecting associated biological ground truth and calibration information at a scale compara-5 Compact representation of echosounder time series ble to that of the acoustic data. For example, extensive net-or optics-based biological taxa identification samples are critical in constructing empirical reference or training datasets for scatterer classification, but the cost and logistical complexity involved in collecting them in parallel with autonomous echosounder sampling is often prohibitive, especially in geographically remote areas or deep water. The validity and applicability of any derived classification rules further hinge on the assumed invariance of the combination of measurement systems and the biological community, both of which can vary significantly depending on time and space. In addition, accurate echosounder calibration-the mapping between the received voltage signal and the actual echo strength in physical units (Demer et al., 2015)-imposes a stringent requirement for adequate application of these methods, as calibration needs to be conducted separately for each frequency under a wide range of environmental conditions (e.g., at depth) under which acoustic data are collected.
Another suite of methods were developed with a goal of capturing specific biological activity patterns observed in the echogram. These methods seek to capture local echogram features through summary statistics derived from a vertical echogram slice or across a few immediately adjacent slices. The features include integrated echo energy (as a proxy for biomass), the vertical spread of animal distribution, timing and speed of animals engaged in vertical migration, and the number and depth of strong scattering layers, etc. (e.g., Cade and Benoit-Bird, 2014;Klevjer et al., 2016;Proud et al., 2015;Urmy et al., 2012;Woillez et al., 2007). While useful in summarizing echogram variations, these statistics are handcrafted and not adaptive, risking loss of information or poor characterization when the pre-determined echogram features mismatch those of the actual data. The summary 6 Compact representation of echosounder time series statistics also lack the capability to automatically extract and describe high-level echogram features since the computation involves only temporally and spatially localized measurements in large datasets. A recent study circumvented some of these difficulties by employing a timeaugmented extension of the empirical orthogonal function method to examine the anomaly of zooplankton diel vertical migration behavior with respect to monthly means (Parra et al., 2019). The anomaly patterns were subsequently linked to factors such as ocean currents, lunar variability and the taxonomic structure of the local biological community.
Rather than focusing on specific organisms or localized echogram features, the overarching goal of our work is to develop a methodology that can automatically build a comprehensive representation of high-level spatio-temporal structures intrinsic to the echosounder data, such as temporally transient appearance of scattering layers or animal aggregations engaged in different vertical movement behaviors. These structures are the most prominent and immediate features human observers attend to in visual analysis of large echograms, and are key to forming and testing hypotheses pertinent to organism responses to environmental disturbance of natural or anthropogenic sources.
In this paper, we show that matrix decompositions are powerful techniques capable of achieving this goal in an unsupervised manner with minimal assumption and prior knowledge from the observed region. We show how different decomposition formulations are suitable for exploiting and extracting different structures in the data: 1) Principal Component Pursuit (PCP) is a robust version of Principal Component Analysis which exploits the latent regularities in long-terms echograms to automatically remove spurious and noisy echo outliers while retaining the high-level spatial and temporal patterns present in the original time 7 Compact representation of echosounder time series series; 2) temporally smooth Nonnegative Matrix Factorization (tsNMF) is a decomposition with nonnegativity and temporal smoothness constraints, which successfully extracts distinct daily echogram patterns with time-varying activation sequences, providing a compact representation of temporal processes embedded in the data. The utility of these methods 8 Compact representation of echosounder time series are demonstrated using a three-frequency echosounder time series spanning 62 days. While no specific constraints are incorporated into the decomposition formulation in this paper to account for multi-frequency echo features, we will discuss preliminary results from a recent work in which we investigated how these features may be exploited through nonnegative tensor decomposition, a multi-way extension of nonnegative matrix factorization as well as avenues for future development (Lee and Staneva, 2019).
Our contribution builds on the classic concept of dimensionality reduction to find structures directly from the data and use them for representation and interpretation. In contrast to previous methods that depend strongly on past experience and observer assumptions, the unsupervised data-driven approach we employ is adaptive to intrinsic properties of the data. Such a capability is pivotal in extracting and synthesizing information from the overwhelming volumes of echosounder data collected from not only unfamiliar ecosystems but also those undergoing rapid changes in the changing climate. The raw data collected by the echosounder are preprocessed using the Python package echopype (Lee et al., 2020) prior to the decomposition analyses. The preparation steps include: 1) parsing and saving the .raw binary data files to self-described machine-readable netCDF files, 2) performing calibration and echo-integration to obtain volume backscattering strength (S v , units: dB re 1 m −1 ) for each echo sample (each ping) (Simmonds and MacLennan, 2007), and 3) obtaining mean volume backscattering strength (MVBS or S v ) over non-overlapping echogram regions with a depth bin size of 5 m and a temporal bin size of 200 sec. Since this echosounder was not formally calibrated during this deployment, we use manufacturer-supplied parameters stored in the data files to calibrate the echo measurements.

A. Acoustic and environmental data
In addition to data from the echosounder, ocean current measurements collected by an acoustic Doppler current profiler (ADCP) from the OOI Oregon Offshore Cabled Benthic Experiment Package (Ocean Observatories Initiative CE04OSBP) are also obtained for joint analysis of the decomposition results.
10 Compact representation of echosounder time series  (Bouwmans et al., 2016;Cichocki et al., 2009). These methods leverage the hidden low-dimensional nature of high-dimensional observations and yield latent factors that are usually more tractable and interpretable compared to the original data. These techniques have found wide applications in many scientific and applied fields, including computer vision, neuroscience, text mining, deep learning, etc. (Chen and Cichocki, 2005;Cichocki et al., 2009;Ding et al., 2008;Lee and Seung, 1999;Oseledets, 2011;Pnevmatikakis et al., 2016). A notable and illustrative application example is the analysis of human faces using Principal Component Analysis (PCA), which produces a series of "eigenfaces" that capture increasingly sophisticated variation of facial features across subjects (Turk and Pentland, 1991).
Using matrix decomposition methods, here we aim to automatically extract high-level spatio-temporal structures embedded in the echogram of each day and use these patterns as descriptors for temporal changes in the time series. The choice of the daily interval is motivated by the profound influence of daylight cycle to life in the ocean and the significant contribution of diel vertical movements of marine animals in the global biogeochemical cycle (Hays, 2003;Lehodey et al., 2010). Below we describe the data restructuring procedure and the decomposition formulations used in this study to remove noisy outliers from data, discover distinct echogram patterns, and summarize variations in long-term echosounder time series.

12
Compact representation of echosounder time series

Restructuring echo data for decomposition
The multi-frequency echosounder time series are restructured to allow two-way (matrix) decomposition in this paper. All echo observations from a single frequency of each day are first "flattened" and concatenated to form a long vector that contains all echogram pixels within the day. The long vectors from all days in the observation period are then combined horizontally along the day dimension to form a two-dimensional matrix (Fig. 3A) for each frequency. The matrices from all frequencies are then stacked vertically to form the final data matrix for decomposition (Fig. 3B). Each single-frequency daily echogram is of dimension is the number of averaged depth bins and N p (=144) is the number of averaged ping bins of the MVBS data within a day. Therefore, the "flattened" long vectors is the number of frequencies and N t (=62) is the number of days.

Outlier removal
A profiler collocated with the echosounder on the same mooring traversed the water column at irregular time intervals during the observation period and introduced strong interference to echoes from natural sources (indicated by arrows in Fig. 2A). The profiler echo interference needs to be removed from the time series before decomposition analysis to avoid their dominance in the optimization updates. To this end, we employ Principal Component Pursuit (PCP), or robust PCA (Candes et al., 2009), which was designed to separate the 13 Compact representation of echosounder time series Here we use PCP as a purely data-driven approach to isolate abrupt patches of strong echoes from the profiler and fine echogram granularity (the sparse component) from highlevel spatio-temporal structures in the echogram (the low-rank component), which are the focus of subsequent analyses. The decomposition was performed on the restructured MVBS matrix described in Sec. II B 2. The PCP decomposition is performed using the Python package RobustPCA (Shun Chi, 2018).
Beyond the scope of this paper, PCP may be useful in removing electric or acoustic interference commonly observed between unsynchronized active acoustic instruments on ships (e.g., ADCP and echosounders), and in isolating echoes from dense schools of pelagic fish that tend to be spatially and temporally localized on the echogram.

Additive representation of temporal processes in a long-term echogram
Nonnegative Matrix Factorization (NMF) is a data decomposition alternative to PCA for which both the latent factors and their activation coefficients are constrained to be nonnegative (Lee and Seung, 1999;Paatero and Tapper, 1994) (we refer this original formulation as the "traditional NMF" hereafter). The nonnegativity constraints make the latent factors (components) purely additive building blocks that can be combined to reconstruct the data. Moreover, the components tend to contain primarily localized features, which underlie a more interpretable "part-based" reconstruction and representation. The advantage is readily demonstrated by the decomposition of human faces into parts associated with eyes, 15 Compact representation of echosounder time series mouths, etc. shown in Lee and Seung (1999), as opposed to the eigenfaces from PCA that each contains variations of different parts of the entire face and can sometimes be challenging to interpret.
In this paper we use a regularized temporally smooth Nonnegative Matrix Factorization (tsNMF) formulation to discover patterns (components) embedded in the low-rank MVBS echogram and use the activation of these patterns to describe temporal variations in the echosounder time series. By performing decomposition along the dimension of day (see and space (depth) that co-vary throughout the time series. These spatially and temporally coupled "daily echogram patterns" are ecologically relevant, as organisms in each group are bound together by their movement behavior and not necessarily taxonomic identity.
tsNMF is a modification of traditional NMF and incorporates in its formulation a temporal smoothness constraint which controls abrupt changes in the pattern activation sequences and emphasize the role of sequential information in the decomposition (see below).
We optimize the following cost function in tsNMF as proposed by Fabregat et al. (2019): where X is the data matrix (the low-rank MVBS echogram L from PCP in this study), W contains the latent factors (patterns), H contains the patterns' activation coefficients, T is 16 Compact representation of echosounder time series the Tikhonov matrix that regularizes the smoothness of H through the parameter η, λ allows tuning for sparsity in W, β W and β H avoid matrices W and H to grow too large, and i, j, k are indices for the matrix elements. · F is the Frobenius norm for matrices and is equal to the square root of the sum of squares of all matrix entries, while · 1 is the l 1 matrix norm defined as In the context of this paper, each column w in W encodes a particular daily echogram pattern and each row h in H captures the activation of the corresponding pattern in the observation period (Fig. 3B). Since MVBS values in the data matrix are negative in absolute physical units, we perform nonnegative decomposition by first subtracting the minimum of the entire matrix and adding the minimal value back for reconstruction. This approach is justified by our primary focus in extracting high-level patterns in the data, rather than physics-based inversion of scatterer numerical density (see Sec. IV for more discussion).

Selection of decomposition parameters
The decomposition parameters, including rank, regularization parameters and initialization of H and W are selected based on empirical analysis of the decomposition outcome.
Finding the intrinsic low dimension of a dataset has long been a challenging task for many unsupervised techniques, including clustering (Mirkin, 2011). Here we use the change of reconstruction error across increasing rank (Frigyesi and Höglund, 2008;Hutchins et al., 2008) and the cophenetic correlation coefficients (Brunet et al., 2004) as guides for rank selection.
The former provides an indication of the rank beyond which no significant non-random structures in the data are captured by the decomposition, and the latter provides a measure for decomposition stability across multiple runs. The smoothness parameter η is selected based on the L-curve concept, which seeks to balance the reconstruction error and the smoothness regularization cost [the first and the second term, respectively, of eq. (1)], which co-vary along an L-shaped curve that contains optimization regimes dominated by either the former or the latter term as the smoothness parameter sweeps across a wide value range (Cheng 18 Compact representation of echosounder time series cluster formed by these two terms over 320 randomly initialized runs.

C. Summarizing variations in a long-term echosounder time series
The activation of daily echogram patterns from matrix decomposition provides a compact representation that is useful for summarizing temporal processes in a long-term echosounder time series. To assess changes in the overall echogram structure, we compute the Euclidean distance between the pattern activation for pairs of days and perform hierarchical agglomerative clustering using linkage calculated using Ward's minimum variance method. The distance matrix and clustering outcome allow convenient identification of the transitions of overall echogram structure throughout the observation period.

19
Compact representation of echosounder time series

III. RESULTS
A. Outlier-free low-rank echogram PCP successfully decomposes the MVBS echogram into a low-rank component that contains the high-level spatio-temporal features in the data (Fig. 2B) and a sparse component consisting of abrupt echo patches and fine-grained details (Fig. S1). The low-rank component series. This is achieved without prior knowledge of the profiler movement schedule, which was in fact irregular during the observation period due to maintenance needs.

B. Daily echogram patterns as building blocks of long time series
The low-rank MVBS echogram from PCP is decomposed by tsNMF (rank=3, η = 500000) into three distinct daily echogram patterns with corresponding activation sequences throughout the observation period (Fig. 4). Three is the lowest rank at which the decomposition produces clearly distinguishable daily echogram patterns and at the same time captures the majority of the structures in the data (for additional details see Fig. S2 and S3).
The three daily echogram patterns collectively represent a repertoire of daily echo energy distribution patterns that jointly reconstruct the low-rank MVBS time series, demonstrating the "parts-based representation" property of nonnegative decomposition (Fig. 4, Sec. II B 4).
Pattern #1  Pattern #3 represents an aggregation of scatterers that appear in the mid to upper water column during daytime, with a trend of decreasing echo strength over increasing frequency that is consistent with fish-like scatterers. While its activation is more difficult to discern than the other two patterns, its inactivation produces noticeable changes on the echogram: the upper water column appears significantly "emptier," especially at 38 kHz, when Pattern #3 activation is low around August 22 and from late September to early October.

23
Compact representation of echosounder time series during the 62-day observation period (Fig. 5A). Along the matrix diagonal, groups of consecutive days that share similar co-activation of echogram patterns are clearly identifiable as reddish "blocks." The off-diagonal entries show similarity between temporally more distant observation days and reveal recurring co-activation of the daily echogram patterns. The echogram structural changes are also illustrated by the multiple transitions between four clusters of pattern co-activation identified by agglomerative clustering based on pairwise distance between days ( Fig. 5B and Fig. S4). Pairing the distance matrix and the cluster transition diagram with ADCP measurements further reveals interesting coincidence at multiple time points between overall echogram structural changes and prominent direction reversal or magnitude variation of ocean currents at the upper 300 m of the water column (indicated by the arrows and dash lines in Fig. 5C). For example, the pattern activation cluster transitions at September 3, 10, 22, and 24 coincide with direction reversal in one or both of the eastward and northward components of the currents. The cluster transitions at August 25 and October 7 additionally coincide with an increase of current velocity in the lower and upper water column, respectively.

IV. DISCUSSION
In this paper we present automatic extraction and compact representation of temporal processes in a long-term echosounder time series via matrix decomposition. We show that unsupervised decomposition techniques are capable of transforming complex echo observation into low-dimensional components that are more tractable and interpretable than the original data. The decomposition methods automatically remove noisy outliers and capture high-

24
Compact representation of echosounder time series level spatio-temporal echogram structures through time-varying linear combination of a small number of ecologically relevant daily echogram patterns. The implication of our results is profound: this approach is entirely data-driven and does not assume knowledge of the type of scatterers or sources of interference present in the data. Such an approach is fundamentally different from previous echo analysis methods that rely heavily on handcrafted rules built from past experience, limited biological ground truth information, and observer assumptions.
The methodology we present here is therefore more flexible and adaptable when applied to data collected from geographical regions or environmental conditions for which scarce prior information exists. This type of data are now prevalent due to the broad availability of echosounders-equipped autonomous ocean observing platforms and as a result of the rapidly changing climate.
Our results show that the parameter-free, automatic PCP decomposition significantly reduces the dimensionality of the MVBS echogram (Fig. 2) and enables the subsequent tsNMF decomposition analysis. By restructuring data according to the daily interval (Fig. 3), we exploit the intrinsic regularity in the data through PCP to remove the echo interference from the water column profiler. Without this crucial step, the profiler echo traces would have dominated the tsNMF optimization updates due to its significant echo strength. In addition, while much fine echogram granularity is decomposed into the sparse component by PCP, the low-rank MVBS echogram retains the high-level spatio-temporal structures intrinsic to the original data that are important for joint analysis of echosounder observation with other ocean environmental parameters.

25
Compact representation of echosounder time series The results also show that a temporally regularized nonnegative decomposition through tsNMF provides a tractable and interpretable dissection of the echosounder time series (Fig. 4). The decomposition discovers daily echogram patterns that each represents a distribution of echo energy across time (hour of day) and space (depth), with co-varying activation throughout the observation period. For biological scatterers, such grouping is therefore based on behavior and does not necessarily correspond to specific taxonomic groups of organisms that have been the focus of many echo classification methods (see Sec. I). This behavioral binding can be ecologically significant, since spatial and temporal co-occurrence of organisms within daily movements suggests trophic functional grouping that mediates energy transfer throughout the water column (Lehodey et al., 2010).

Modeling frequency domain information
Interestingly, while frequency domain information is not explicitly modeled in tsNMF, the constraint on temporal smoothness of pattern activation encourages biological (and therefore spectral) coherence within each pattern, since movements of animal aggregations are continuous in time and space. Recall that spectral characteristics of echoes vary significantly depending on the size, shape, orientation, and material properties of the scatterer. In a preliminary study we have examined the effects of explicitly modeling frequency domain information via tensor decomposition, which extends matrix decomposition to a multi-way component analysis (Lee and Staneva, 2019). While the nonnegative tensor decomposition yielded three daily echogram patterns that resemble the general distribution of echo energy in Pattern #1 and #2 discovered by tsNMF here, the spectral dependency within each pat-

26
Compact representation of echosounder time series tern is different. Further investigation revealed that the reconstruction error is substantially higher for tensor decomposition, rendering its results less reliable and less interpretable. This is likely due to the combination of 1) the restrictive form of Kruskal tensor decomposition we used, which requires coherent variation across frequency for all pixels within a single pattern but at the same time does not allow interactions across patterns in the reconstruction (Cichocki et al., 2009), and 2) the small number of frequency features (N f = 3) in the dataset, which limits the total number of patterns (rank) allowed in the decomposition. Future investigations using data with richer spectral features, such as those collected by broadband echosounders, would be more appropriate for evaluating the utility of tensor decomposition in modeling frequency domain information in the decomposition.

Relation to environmental data
The pattern activation sequences from tsNMF and the derived distance matrix summarize both similarities and changes in the echosounder time series, providing an explicit avenue for joint analysis of acoustic observation and ocean environmental data. By aligning ocean currents measured by ADCP with the distance matrix and activation clusters, a clear pattern emerges that associates changes in the echogram structure to changes of the direction or magnitude of ocean currents at multiple time points throughout the observation period (indicated by arrows in Fig. 5C). These associations are ecologically plausible, since midtrophic animals can be advected by strong currents, and their vertical movement behaviours may also be modulated by changes in the ambient environment. While detailed interpretation of these observations is beyond the scope of this paper, we note that strong currents (up to 27 Compact representation of echosounder time series 0.2 m/s) could be biologically significant in this region (Keister et al., 2009;Wu et al., 2014).
Using an example dataset consisting of 62 days of observation, we illustrate the synoptic summarizing power of unsupervised decomposition that will likely play a crucial role in extracting information from much longer time series that are difficult to analyze manually.

Comparison with PCA
Our results demonstrate the advantage of interpretability from imposing nonnegative constraints in decomposition analysis of echosounder data. In comparison with PCA, the purely additive linear combination of low-rank patterns makes it straightforward to attribute changes in the long-term echogram to temporally varying contributions from different daily echogram patterns. For example, while the appearance of the first two PCA components resembles the first two tsNMF patterns, their contributions to the observed MVBS echogram are much more difficult to interpret due to the intermingled positive and negative entries in both the patterns and the activation sequences (Fig. S5). PCA could however useful for analyzing anomaly with respect to a mean pattern, as has been shown in Parra et al. (2019).

Rank and parameter selection
The selection of hyperparameters in the decomposition, such as the rank and the smoothness parameter, is a topic that requires attention in future application of similar methods.
In this work we use intrinsic properties of the data to select these hyperparameters: we select the minimum rank that yields distinctive daily echogram patterns with moderate re-  . This comparison provides a guide for selecting the rank beyond which no significant data structures are captured by the decomposition (Frigyesi and Höglund, 2008;Hutchins et al., 2008).
The median MSE of all 320 randomly initialized runs are plotted. Each curve is normalized by subtracting the medium MSE at rank=1. (B) Variation of the cophenetic correlation coefficient with increasing rank. This coefficient provides a measure for decomposition stability (Brunet et al., 2004) and is calculated by grouping days according to the daily echogram pattern with the highest activation strength (the "connectivity") and evaluating the consensus of the connectivity over repeated runs. Here the consensus is computed as a weighted average based on the MSE following the modification proposed by Frigyesi and Höglund (2008). Note that the computation of connectivity requires assigning the membership of a particular observation day to only one echogram pattern, which conflicts with the fact that it is the joint activation of multiple patterns reconstructs the data. Therefore, we base our rank selection primarily on MSE variation and use the variation of cophenetic correlation coefficient only as a reference.

41
Compact representation of echosounder time series is changed to a bi-directional diverging scheme to highlight the fact that PCA components contain both positive and negative entries, instead of the solely nonnegative entries from tsNMF. While the first two PCA components resemble the first two from tsNMF (Fig. 4), their contributions to the observed MVBS echogram are much more difficult to interpret due to the intermingled positive and negative entries in both the patterns and the activation sequences.