Machine learning technique using the signature method for automated quality control of the Argo profiles

TL;DR

This paper introduces a machine learning approach using the signature method to automate quality control of Argo ocean profiles, offering a simple, objective alternative to rule-based methods by leveraging the mathematical properties of path signatures.

Contribution

The study applies the signature method to Argo profile data for supervised learning, demonstrating its effectiveness in automatic quality control without subjective rule-based criteria.

Findings

Achieved accurate quality control with cross-validation.

Provided a mathematically grounded, objective method.

Outperformed traditional rule-based approaches.

Abstract

A profile from the Argo ocean observation array is a sequence of three-dimensional vectors composed of pressure, salinity, and temperature, appearing as a continuous curve in three-dimensional space. The shape of this curve is faithfully represented by a path signature, which is a collection of all the iterated integrals. Moreover, the product of two terms of the signature of a path can be expressed as the sum of higher-order terms. Thanks to this algebraic property, a nonlinear function of profile shape can always be represented by a weighted linear combination of the iterated integrals, which enables machine learning of a complicated function of the profile shape. In this study, we performed supervised learning for existing Argo data with quality control flags by using the signature method, and demonstrated the estimation performance by cross-validation. Unlike rule-based approaches, which require several complicated and possibly subjective rules, this method is simple and objective in nature because it relies only on past knowledge regarding the shape of profiles. This technique should be critical to realizing automatic quality control for Argo profile data.

Tables

Table 1. Table 1: Confusion matrix with cutoff y c subscript 𝑦 𝑐 y_{c}

	$0$	$1$
	True-positive $N_{\text{TP}}$	False-positive $N_{\text{FP}}$
$y_c\le \tilde{y}$	False-negative $N_{\text{FN}}$	True-negative$N_{\text{TN}}$

Equations

$${\mathrm{X}}^{(i_1\cdots i_n)}$$

$${\mathbf{X}}_n$$

$${\mathcal{S}}^n(X)$$

$${\mathcal{S}}^2(X)$$

$${\left(X\ast Y\right)}_{\tau }=\{\begin{array}{ll}{X}_{\tau }& \text{if }\tau \in [0,s]\\ X_s+Y_{\tau }-Y_s& \text{if }\tau \in [s,t].\end{array}$$

$${\mathcal{S}}^n(X)\otimes {\mathcal{S}}^n(Y)$$

$${\mathbf{Z}}_h$$

$${\mathrm{Z}}^{(i_1\cdots i_h)}$$

$${\mathcal{S}}^2(X)\otimes {\mathcal{S}}^2(Y)$$

$${\mathcal{S}}^n(X\ast Y)={\mathcal{S}}^n(X)\otimes {\mathcal{S}}^n(Y)$$

$${\mathrm{X}}^{(P)}$$

$${\mathrm{X}}^{(PP)}$$

$${\mathrm{X}}^{(i)}$$

$${\mathrm{X}}^{(ijk)}$$

$${\mathcal{S}}^n(X_{\text{seg}})$$

$${\mathcal{S}}^n(X_{s,t})$$

$$X$$

$$\displaystyle\begin{bmatrix}\color[rgb]{0,0,1}{P_{0}}\\ \color[rgb]{0,0,1}{S_{0}}\\ \color[rgb]{0,0,1}{T_{0}}\\ \color[rgb]{0,0,1}{P_{0}}\\ \color[rgb]{0,0,1}{S_{0}}\\ \color[rgb]{0,0,1}{T_{0}}\end{bmatrix},$$

$$y=\sum _I{w}^I{\mathrm{X}}^I+ϵ,$$

$${\mathrm{X}}^A{\mathrm{X}}^B$$

$$J(w)$$

$$\tilde{y}(m)$$

$$\mathrm{X}(m{)}^I$$

$$\tilde{y}(m)$$

$$U=\left[\begin{array}{ccc}{U}_1^1& \cdots & U_1^L\\ ⋮& & ⋮\\ U_{N_{\text{pc}}}^1& \cdots & U_{N_{\text{pc}}}^L\end{array}\right]$$

$$r(m)$$

$$d{Y}_{\tau }^i$$

$$Y_t^i$$

$$\varphi \left(\{X_{\tau }{\}}_{0\le \tau \le t}\right)$$

$$f\frac{\partial u}{\partial P}=-{\frac{\partial }{\partial y}({\rho }^{\bullet -1})|}_{P=\text{const.}},f\frac{\partial v}{\partial P}={\frac{\partial }{\partial x}({\rho }^{\bullet -1})|}_{P=\text{const.}},$$

$$f{\int }_{x^{\prime }=x_0}^{x_1}\frac{\partial v}{\partial P}d{x}^{\prime }$$

Full text

\draftfalse\journalname

Earth and Space Science

Research and Development Center for Global Change, JAMSTEC, Yokosuka, Japan

\correspondingauthor

Nozomi Sugiuransugiura@jamstec.go.jp

{keypoints}

Machine learning for the quality control flags of Argo profiles was performed.

By converting each profile sequence of temperature, salinity, and pressure into its signature, classification was performed efficiently.

The signature is regarded as a fundamental object that represents a data sequence.

1 Introduction

Argo is an international effort collecting high-quality temperature and salinity profiles, typically from the upper 2000 m of the global ocean [Gould \BOthers. (\APACyear2004)]. The data come from battery-powered autonomous floats that drift mostly at a depth where they are stabilized at a constant pressure level. At typically 10-day intervals, the floats rise to the surface for approximately 6 h while measuring temperature and salinity. On surfacing, the satellites position the floats and receive the transmitted data. Now, the array of over 3000 floats provides 100,000 temperature/salinity profiles annually distributed over the global oceans at an average 3-degree spacing. The quality control (QC) of the massive Argo profile data [ARGO (\APACyear2019)] must be systematic to keep the quality of the observational data homogeneous and to utilize human resources efficiently. In addition, accurately quantifying the relationship between the profile shape and the effect it has on oceanic processes is essential for understanding the ocean state through the profile observation. Conventionally, significant time and effort are spent assigning the quality control flag to each Argo profile.

Regarding attempts for advanced automatic QC procedures on oceanographic profiles, some studies have been applied to the Argo CTD profile because of the huge amount of data accumulated for 200 million profiles over 20 years. For example, \shortciteAUdaya2013 provided a semi-automatic QC procedure using objective mapping to remove anomalous values from the profiles. \shortciteAUDAYABHASKAR2017469 demonstrated automatic QC by defining the convex fulls from the climatological dataset. Meanwhile, \citeAOno2015 attempted to apply a machine learning method to the delayed-mode QC of Argo profiles towards a possible automatic QC system for an Argo data stream. Similarly, an integrated Argo data flow using machine learning was introduced to be an automated system with an improved QC ability (presented by \shortciteAMaze2017 in the report of the 18th Argo Data Management Meeting). Thus, QC procedures for oceanographic data have been gradually improved by many researchers using advanced tools or methods.

The discrimination procedures involved in the automation of the quality control have been performed mainly in a rule-based manner [<]e.g.¿[]7344896,7838290,7849862. As an alternative and more flexible approach, this study attempted to automate the process via supervised learning of the human judgment process. In doing so, it is essential to quantify the profile shape so that the function that yields the quality control flag can be expressed as a linear combination of the numerical values that represent the profile shape. The machine learning thereby reduces to a linear optimization problem that can be easily solved. The key tool that enables this quantification is the signature, which is the set of all iterated integrals [Chevyrev \BBA Kormilitzin (\APACyear2016), Levin \BOthers. (\APACyear2013)], proposed in the theory of rough path by \citeAlyons2007differential.

In this research, we propose a procedure of first converting the vector sequence of each Argo profile into a sequence of real numbers that represents its shape and then expressing a nonlinear function of the shape in the form of a linear combination of these numbers; this conversion facilitates machine learning of the nonlinear function. A machine learning experiment regarding the function was performed and applied to automatic assignment of quality control flags to the profiles.

2 Theoretical background

The central concept in this study is the signature, proposed in the theory of rough path by \citeAlyons2007differential. In what follows, we briefly introduce the concept of signature and the notation used in this paper. For more details, refer to \citeAchevyrev2016primer,

As perceived from a re-examination of controlled differential equations (refer A), characteristics of a data sequence can be represented by the signature, which comprises the iterated integrals. Note, in this paper, the subscript notation $X_{\tau}$ is used to denote dependence on the parameter $\tau\in[0,t]$; $A^{\bullet n}$ and $A^{\otimes n}$ denote the $n$-th power and $n$-times tensor product, respectively, but otherwise, a superscript denotes a component.

Suppose we have a sequence of $d$-dimensional vectors $X_{u}~{}(0\leq u\leq t)$. Let the time order be $0<t_{1}<\cdots<t_{n}<t$. We define the iterated integral for indices $i_{1},\cdots i_{n}=1,\cdots,d$ as

[TABLE]

where we should be careful about the difference between the font for a sequence of vectors $X^{k}_{t_{1}}$ and the one for an iterated integral $\mathrm{X}^{(k_{1}k_{2}\cdots k_{n})}$. By treating all the index values together, we obtain a tensor of order $n$:

[TABLE]

and $\mathbf{X}_{0}$ is constant $1$. Moreover, by putting together the iterated integrals for all combinations of the indices, we obtain the signature up to degree $n$:

[TABLE]

which has $(d^{n+1}-1)/(d-1)$ components. For instance, the signature up to degree $2$ for a $2$-dimensional sequence is

[TABLE]

where $X_{0,t}^{i}\stackrel{{\scriptstyle\text{def}}}{{=}}X_{t}^{i}-X_{0}^{i}$, $\bullet 2$ denotes the second power, and $\mathrm{X}^{()}=\mathbf{X}_{0}=1$. Note that the order of integrands matters in $\int_{0}^{t}X_{0,u}^{1}dX_{u}^{2}$ and $\int_{0}^{t}X_{0,u}^{2}dX_{0,u}^{1}$. In general, it is important for the signature to encode the order in which each component changes along the path.

Suppose we have two paths, $X=\{X_{\tau}\}_{0\leq\tau\leq s}$ and $Y=\{Y_{\tau}\}_{s\leq\tau\leq t}$. Their concatenation is the path defined by

[TABLE]

On the other hand, regarding their signatures, $\mathcal{S}^{n}(X)=(\mathbf{X}_{0},\mathbf{X}_{1},\cdots,\mathbf{X}_{n})$ and $\mathcal{S}^{n}(Y)=(\mathbf{Y}_{0},\mathbf{Y}_{1},\cdots,\mathbf{Y}_{n})$, we can define the product as

[TABLE]

whose components are

[TABLE]

For instance, the product of the signatures, up to degree $2$, for the $2$-dimensional sequence is

[TABLE]

In this manner, the set of signatures has a group structure in the free tensor algebra with respect to the product $\otimes$. Furthermore, Chen’s identity [Chen (\APACyear1958)]:

[TABLE]

is satisfied, which defines a homomorphism from path space with concatenation (5) to signature space with group operation (6).

In the context of geophysics, we can show that some diagnoses for oceanographic conditions are written in terms of iterated integrals. Consider a vertical sequence of vector $(P_{\tau},S_{\tau},T_{\tau})$ (pressure, salinity, and temperature) in the ocean.

The first-order iterated integrals are

[TABLE]

which are profile depth, sea surface salinity, and sea surface temperature, respectively. 2. 2.

The second-order iterated integrals include

[TABLE]

which represent the square of profile depth, total salinity content, and total heat content, respectively.

We find another example in B.

Note that $P$ is treated equally to $S,T$ in the above because the seemingly redundant parameter $\tau$ is essential to ensure that the path has no self intersection and the signature is invariant under the reparameterization of $\tau$. If one parameterizes $T,S$ with $P$, the path would be drawn on a two-dimensional $T,S$-surface, which loses considerable information on the shape of the sequence.

3 Method

The data used in this research were observed by the global array of Argo floats [ARGO (\APACyear2019)], each of which floats and sinks from the sea surface to a depth of approximately $2000\mathrm{m}$.

Because the shape of a vector sequence $(P_{\tau},S_{\tau},T_{\tau})$ is only perceived in a certain reference frame, it is convenient to make the original quantities dimensionless; $P$ in $\mathrm{dbar}$, $S$ in $\mathrm{psu}$, and $T$ in °C into $\widehat{P}=P/2000$, $\widehat{S}=S/2$, and $\widehat{T}=T/20$, where divisor $(2000,2,20)$ is chosen as a typical scale of the components. For simplicity, henceforth, we omit the hat symbol for the component. Figures 1 and 2 show examples of the vertical profiles of temperature, salinity, and pressure, along with the corresponding iterated integrals. By virtue of quality control procedures with manual judgment, the quality control flags are already assigned to all of the data.

Here, we describe the basic concept of the signature method and how to apply it to Argo profiles. We also explain how to construct a procedure for supervised learning using the signature and how to verify the results.

3.1 Representing the Argo profile shape by signature

3.1.1 Computation of signature

Suppose we have $d$-dimensional profile data $\{X_{\tau}\}_{0\leq\tau\leq t}$ that can be seen as a line graph connecting points $X_{0}=X_{u_{1}},X_{u_{2}},\cdots,X_{u_{L}}=X_{t}$; then, we can compute its iterated integrals as follows:

For line segment $X_{\text{seg}}\stackrel{{\scriptstyle\text{def}}}{{=}}\{X^{i}_{u}+X^{i}_{u,u^{\prime}}\tau\}_{0\leq\tau\leq 1}^{i=1,2,\cdots,d}$, which has starting point $X^{i}_{u}$ and slope $X^{i}_{u,u^{\prime}}$, the iterated integrals are calculated as

[TABLE]

and the [math]-th iterated integral is constant $1$. In this case, the signature (up to degree $n$) is nothing but a commutative exponential function for the vector $X_{u,u^{\prime}}$:

[TABLE]

where $\mathbf{e}_{i_{k}}$ is the $i_{k}$-th unit vector. 2. 2.

Let the time order be $s\leq u\leq t$. By concatenating a path $X_{s,u}$ from time $s$ to $u$ with a path $X_{u,t}$ from time $u$ to $t$, we obtain a path $X_{s,t}$ from time $s$ to $t$, whose signature is the product of the signatures:

[TABLE]

which is due to Chen’s identity (10). 3. 3.

By concatenating the paths successively using Eq. (13), we can compute the signature for the whole line graph.

The numerical computation of the signature in this study is performed by using Python library Esig [Kormilitzin (\APACyear2017)].

3.1.2 Lead-lag transformation

Suppose we have a sequence of $d_{0}$-dimensional ($d_{0}=3$) vectors with length $L+1$:

[TABLE]

To more precisely grasp the shape of the line graph, we perform a lead-lag transformation [Chevyrev \BBA Kormilitzin (\APACyear2016)], which defines a sequence of $d(=2d_{0})$-dimensional vectors with length $L\cdot d+1$:

[TABLE]

The transition rule for the lead-lag transformation is as follows:

Take two copies of $X_{0}$ and use it as the initial condition. 2. 2.

Update only $1$ component among $d$ components at once. 3. 3.

Use the previous value instead if the present value is missing.

3.2 Machine learning procedure for quality control process

Suppose we have a set of profile data $X(m)\stackrel{{\scriptstyle\text{def}}}{{=}}\{X_{\tau}(m)\}_{0\leq\tau\leq t}$ for $m=1,2,\cdots,M$, whose signature is denoted as $\mathrm{X}(m)\stackrel{{\scriptstyle\text{def}}}{{=}}\mathcal{S}(X(m))$. Let us consider the problem of assigning the discriminant values to each profile depending on whether a profile matches the quality standard.

We first make a model for the rule of quality control as a functional form; that is, a linear combination of the iterated integrals $\mathrm{X}^{I}$ for all combinations of indices $I=(),(i_{1}),(i_{1}i_{2}),\cdots,(i_{1}\cdots i_{n})$ yields the discriminant value.

[TABLE]

where $\epsilon$ is the error. Since each index in $I$ runs over $1,2,\cdots,d$, the sequence of iterated integrals in Eq. (14) has $\sum_{j=0}^{6}d^{\bullet j}=(d^{\bullet n+1}-1)/(d-1)$ terms. Note that $\mathrm{X}^{()}$ represents the constant $1$. Such a representation is possible because its nonlinearity is unraveled thanks to the property of shuffle product; for a fixed path $X$, the product of iterated integrals for indices $A$ and $B$ is expressed by the iterated integral with respect to the shuffle product $A\mathbin{\mathchoice{\kern 0.0pt\hbox to0.0pt{\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hss}\raise 0.0pt\hbox to0.0pt{\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hfill\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hfill\vrule height=0.0pt,depth=0.0pt,width=0.0pt}\kern 0.0pt}{\kern 0.0pt\hbox to0.0pt{\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hss}\raise 0.0pt\hbox to0.0pt{\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hfill\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hfill\vrule height=0.0pt,depth=0.0pt,width=0.0pt}\kern 0.0pt}{\kern 0.0pt\hbox to0.0pt{\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hss}\raise 0.0pt\hbox to0.0pt{\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hfill\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hfill\vrule height=0.0pt,depth=0.0pt,width=0.0pt}\kern 0.0pt}{\kern 0.0pt\hbox to0.0pt{\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hss}\raise 0.0pt\hbox to0.0pt{\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hfill\vrule height=0.0pt,depth=0.0pt,width=0.0pt\hfill\vrule height=0.0pt,depth=0.0pt,width=0.0pt}\kern 0.0pt}}B$:

[TABLE]

For example, $\mathrm{X}^{(aa)}\mathrm{X}^{(b)}=\mathrm{X}^{(baa)}+\mathrm{X}^{(aba)}+\mathrm{X}^{(aab)}$. This means that a product of iterated integrals is always reduced to the sum of higher-order iterated integrals. Moreover, by virtue of the Stone–Weierstrass theorem, any nonlinear function of the shape of a path can be represented as a linear combination of the iterated integrals. 2. 2.

Suppose we have pairs $(X(m),y(m))$, where each $X(m)$ is a profile sequence, and $y(m)=0,1$ is the discrimination value, which is already given to each sample $m=1,2,\cdots,M$ as training data. Learning these data is simply deriving the weights $w^{I}$ that minimize an $L_{1}$-regularized cost function:

[TABLE]

Because the terms in $\sum_{I}|w^{I}|$ are not quadratic but linear, they have the effect of selecting significant terms under the summation over the set labeled by $I$. This is the notion of least absolute shrinkage and selection operator (LASSO)\shortcite10.2307/2346178 , which can help prevent overfitting. The larger the value of $\alpha$ is, the smaller the number of selected terms with $w^{I}\neq 0$ is. To set an appropriate number of terms, several values of $\alpha$ will be tested. 3. 3.

Using the coefficients $w$ derived in (16), and substituting into Eq (14) the iterated integrals for a profile not used for training, we obtain $\tilde{y}$, an estimate for $y$, as follows.

[TABLE]

where $w^{I}$’s are estimated from the minimization of cost (16). 4. 4.

The minimization problem is efficiently solved by the coordinate descent (CD) method [Friedman \BOthers. (\APACyear2007)].

For the $L_{1}$-regularization term to apply evenly, each iterated integral $\mathrm{X}^{I}$ is preprocessed by subtracting the ensemble mean $\mu^{I}_{\text{train}}$ of the training ensemble and dividing by the standard deviation $\sigma^{I}_{\text{train}}$ of the training ensemble:

[TABLE]

The same operation is performed for the iterated integrals in cross-validation. The minimization problem is solved by using the Python library scikit-learn [Pedregosa \BOthers. (\APACyear2011)].

3.3 Assessment of learning results

The performance of the binary classifier can be quantitatively assessed by visualizing it with the receiver operating characteristic (ROC) curve [Egan (\APACyear1975)]. We refer to the profiles that pass the quality criterion as negative (normal) $y=1$, and the others as positive (bad) $y=0$. By shifting the cutoff value $y_{c}$, one can count the number of positive ones with $\tilde{y}<y_{c}$ and that of negative ones with $y_{c}\leq\tilde{y}$. Then, the samples fall into the four categories in Table 1.

The true-positive rate is defined as $N_{\text{TP}}/(N_{\text{TP}}+N_{\text{FN}})$, and the false-positive rate as $N_{\text{FP}}/(N_{\text{FP}}+N_{\text{TN}})$. The ROC curve is the two-dimensional plot of false-positive rate versus true-positive rate, by changing the cutoff $y_{c}$. It has better performance if the trajectory approaches the upper left corner. Therefore, the area under the ROC curve indicates the performance.

Note that, to improve the readability of the histograms, we use a modified estimation value:

[TABLE]

where we apply transformation $y\mapsto 1-|1-y|$ so that $\tilde{y}(m)\leq 1$.

3.4 Experiment using PCA

Alternatively, principal component analysis (PCA) [<]e.g.,¿thomson2014 can be applied to represent the normal profiles. In that case, the experiment for estimating the quality control flag is performed as follows.

We apply the same nondimensionalization to the $T$- and $S$-sequences as the signature method, and then perform nearest-neighbor interpolation at points $2000\widehat{P}=5,15,\cdots,1995,$ which are placed every $10\mathrm{dbar}$. Accordingly, the sequences are transformed into a sequence $X=(X_{1},\cdots,X_{L})^{T}$ with $L=400$. 2. 2.

Let $\mathcal{F}\stackrel{{\scriptstyle\text{def}}}{{=}}\left\{X(m)|m=1,2,\cdots,M\right\}$ be the set of all profiles. We randomly choose the training ensemble $\mathcal{F}_{\text{tr}}\subset\mathcal{F}$, which comprises negative (normal) samples $\mathcal{F}_{\text{tr},n}\subset\mathcal{F}_{\text{tr}}$, and positive (bad) samples $\mathcal{F}_{\text{tr},p}\subset\mathcal{F}_{\text{tr}}$. 3. 3.

Training is performed by computing the principal components (PCs) for negative training samples $\mathcal{F}_{\text{tr},n}$. Let

[TABLE]

be the truncated PCs, and $\overline{X}=(\overline{X}_{1},\cdots,\overline{X}_{L})^{T}$ be the ensemble mean for the training ensemble $\mathcal{F}_{\text{tr},n}$. 4. 4.

For the $m$-th profile $X(m)\in\mathcal{F}_{\text{tr}}$ (or $\in\mathcal{F}\setminus\mathcal{F}_{\text{tr}}$ for cross-validation), the mean square residual for representing it by the first $N_{\text{pc}}$-PCs is computed as

[TABLE]

where $I_{L}$ is an $L$-dimensional identity matrix. The estimated $\tilde{y}$ is thereby defined as $\tilde{y}(m)\stackrel{{\scriptstyle\text{def}}}{{=}}-r(m).$ 5. 5.

For a fixed threshold value $y_{c}$, we assign negative to the $m$-th profile if $\tilde{y}(m)>y_{c}$ and positive otherwise. The ROC curve is drawn by plotting false-positive rates versus true-positive rates for various $y_{c}$.”

4 Results and Discussion

We used a dataset observed at the location shown in Fig. 3. Each profile is assigned a delayed-mode QC flag by Japan Agency for Marine-Earth Science and Technology (JAMSTEC). We treated profiles with depth widths (the difference between the minimum and maximum depths) of more than $1000\mathrm{m}$, and each profile had approximately $L\sim 100$ observation points. The number of profiles was $M=8.2\times 10^{4}$, and the training data were randomly chosen from these profiles. After applying the lead-lag transformation, each profile was converted into the signature up to order $n=6$.

An overview of the machine learning results is shown by the histogram of estimated values $\tilde{y}$ for normal samples ($y=1$), and the histogram for bad samples ($y=0$).

Figures 4 and 5 show the histograms when $40\%$ of the data are used for training and the remaining $60\%$ are used for cross-validation. We can see that learning is properly performed because there is little difference between the identification of training data and the cross-validation. In particular, this approach never misidentifies negative (normal) profiles if the appropriate cutoff $y_{c}$ is used, but it may accept positive (bad) profiles with a probability $0.6$ when $y_{c}=0.5$. This property is also reflected in the tendency of the ROC curve (Fig. 6) to be almost tangent to the horizontal axis when $x$ is small, but not tangent to $y=1$ when $y$ is large. The histogram for positive samples has two clear peaks, which suggests that the ambiguity is not caused by the judgment by the machine learning, but by the fact that the original quality control flag had a criterion that cannot be decided only by the shape. For example, the original quality control, encoded in $y$, may have a criterion about deviation from climatological variation, which cannot always be detected from profile shape. Moreover, the original quality control is partly done through visual checking, for which the criteria can fluctuate between checks. Obviously, both are not represented by the signature, which is static and shape-oriented.

Figures 7 and 8 shows the histograms when $2.5\%$ of the data are used for training and the remainder is used for cross-validation. In this case, there is a clear tendency of over-learning, which indicates that the number of learning samples, $2.5\%$, is not sufficient.

The performance of a method can be measured from the area under the ROC curves (AUC). Comparing that for the experiments with various ratios of learning samples, we notice that over-learning occurs when the ratio is less than $20\%$ (Fig. 9).

We also compared the results of the experiments with various weights $\alpha$ of the regularization term by the AUC. The number of terms under the summation over the set labeled by $I$ in Eq. (16) is dependent on $\alpha$. Therefore, if we increase the degrees of freedom of the coefficients $w$ by using a smaller $\alpha$, the performance of the reproduction capability increases. However, the estimation capability begins to saturate at approximately $6700$ degrees of freedom (Fig. 10), where $\alpha=10^{-5}$ is used. At that point, the complexity seems to become appropriate.

To confirm the efficacy of lead-lag transformation, we performed a similar experiment as in Figs. 4 and 5 except without lead-lag transformation. We set $\alpha=10^{-5}$ and the proportion of training data to $0.4$. Figure 11 shows the ROC curves for the experiment. The curves for the case with lead-lag are on the upper-left of those for the case without lead-lag, which indicates that the lead-lag transformation helps improve the estimation of the quality control flag.

As a reference case using another representation of the shape, we performed PCA experiments with $N_{\text{pc}}=50$ and $100$PCs. The proportion of training data is set to $0.4$, the same as for the signature case. Figure 12 depicts the ROC curves for the PCA experiments. Although the PCA method also exhibits a considerable skill, the curves stay to the lower-right of those for the signature method, which indicates that the signature method is more effective than the PCA method in estimating the quality control flag. Meanwhile, the computational cost for the signature method is not significantly higher than that for the PCA method, because the former only additionally requires converting each data sequence into the truncated signature, whose calculation load is low.

Further, comparison with the data from the ARGO intercomparison project is performed. The performances of the real-time assignment of QC flags \shortcitewong2020argo by several institutes are shown in \shortciteAwedd2015argo, when the corresponding assignments by the delayed-mode QC are regarded as the ground truth. Because the sets of profile data differ from those in our case, a direct comparison is not strictly relevant but will still serve as a measure of the performance. Figure 13 shows the false-positive vs true-positive rates for those samples in comparison to the signature case. Apart from the the pressure data, all the points for real-time QC data lie on the bottom-right side of our ROC curve. This suggests that the signature method may assign the QC flags more efficiently than the real-time QC procedure does, provided that the past assignment results are ready for use. Another advantage of the signature method is that it assigns the flags consistently to all the components $(P,S,T)$ with higher reliability.

Overall, we found that machine learning using the signature method can learn the existing quality control flags of Argo profiles and automatically assign the flag to new profiles, but it sometimes overlooks bad samples because of the ambiguity inherent in the original quality control flag. Comparative study shows the signature method has a higher performance for estimating the flag than other conventional methods, including the one with the PCA representation and the operational assignments of real-time QC.

5 Conclusions

In this research, we first demonstrated that the shape of a profile from the Argo ocean observing array can be represented by the iterated integrals. Then, we constructed a model for the function that assigns a quality control flag to the shape of a profile, which is expressed as a weighted sum of the iterated integrals.

We performed supervised learning for the weights using the existing quality control flags for training data, and demonstrated via cross-validation that it has good performance in estimating flags for unknown data.

A comparative experiment using the PCA method showed that the signature method, in combination with lead-lag transformation, outperforms the PCA method in estimating the quality control flag. This suggests the superiority of the signature method compared to the conventional machine-learning technique.

This algorithm can potentially enable automatic assignment of quality control flags to new Argo data. The significance of the algorithm is that it objectively and automatically assigns the quality control flag only on the basis of past knowledge about the quality of data without imposing any ad hoc rules. Hence, it should enable more objective and efficient quality control compared to traditional manual methods or rule-based machine learning.

The signature method is quite effective for expressing the shape of an Argo profile and its nonlinear function quantitatively. The rationale for this advantage is that a nonlinear and complicated function of assigning quality control flags can be transformed into a linear combination of the iterated integrals through algebraic transformation (shuffle product) without introducing any errors. This is superior to conventional multivariate regression models, which approximately regard nonlinear dependencies as linear ones. Along this line, we can express, as a function of signature, not only quality control flags but also any oceanic phenomena.

One application of the signature method is assimilation of the signature of observational data into a general ocean circulation model. For example, we can convert a vertical sequence of observational data and that of model data into iterated integrals. We then construct a cost function that compares the signatures for model and observation, rather than directly comparing the state vectors composed of temperature and salinity at each depth. Although a cost term is for a single horizontal and temporal point, data assimilation, in particular the four-dimensional variational method, can combine the effects from multiple terms via model integration and adjoint integration. By doing so, we gain the advantage that the projection of a vertical profile onto any ocean phenomena attains a linear form, which will result in efficient data assimilation. This is expected because many diagnoses for oceanic conditions are written in terms of iterated integrals, as illustrated in sec. 2 and B.

Appendix A Picard iteration

To understand the notion of signature, consider how the theory of rough path treats a data sequence acting on a system. Suppose we have a system of ordinary differential equations with respect to $Y_{\tau}$ forced by a path $X_{\tau}$:

[TABLE]

where $Y^{j}_{\tau}$ is the $j$-the component of vector $Y_{\tau}$, and $F^{i}_{jk}$ is the $i,j,k$-th component of $3$-dimensional tensor $F$.

Performing the Picard iteration yields a solution:

[TABLE]

where $\mathrm{X}_{n}^{(k_{1}k_{2}\cdots k_{n})}\stackrel{{\scriptstyle\text{def}}}{{=}}\int_{0<\tau_{1}<\cdots<\tau_{n}<t}dX_{\tau_{1}}^{k_{1}}dX_{\tau_{2}}^{k_{2}}\cdots dX_{\tau_{n}}^{k_{n}}$ is a component of the $n$-th iterated integral (2). By omitting the indices, we can simply write the solution as $Y_{t}=\left[\sum_{n=0}^{\infty}F^{\otimes n}\mathbf{X}_{n}\right]Y_{0}.$ Notice that the convergence of the series is guaranteed because the magnitude of each iterated integral is uniformly bounded: $|\mathrm{X}_{n}^{(k_{1}k_{2}\cdots k_{n})}|<\frac{L^{\bullet n}}{n!}$, where $L$ is the path length. This form of solution suggests that the action of $X$ on $Y$ can be well summarized by the iterated integrals, and an approximate solution is reproduced by a truncated series of iterated integrals $(\mathbf{X}_{0},\mathbf{X}_{1},\cdots,\mathbf{X}_{n})$, which is called a truncated signature up to order $n$. The point is that the effect of a forcing on a system is asymptotically approximated by the truncated path signature but not by the partial sequence of state vectors.

It has been proven that a path that never crosses itself, like in the case of Argo profiles, is completely determined by its signature [Hambly \BBA Lyons (\APACyear2010)]. A function of a path, say $\phi$, can thus be regarded as that of its signature and compactly approximated by that of a truncated signature:

[TABLE]

A further advantage of such treatment is that the function $f$ can always be expressed as a linear combination of iterated integrals, owing to the shuffle-product property, which is explained later.

Appendix B Thermal wind flow in terms of iterated integrals

As an example of higher-order iterated integrals, we show here that thermal wind flow can be written with iterated integrals.

The thermal wind relation is written in vertical $P$-coordinates as

[TABLE]

where $f$ is the Coriolis parameter, $u,v$ are velocity, $\rho$ is density, and $x,y$ are the longitudinal and latitudinal coordinates, respectively. For a fixed latitude $y$, by performing integrations along the $x$ direction and then the $P$ direction, we obtain an estimate for the meridional velocity as

[TABLE]

where we set $v(x^{\prime},P_{0})=0$ as the layer of no motion. Integrating again along the $P$ direction, we obtain the meridional flow rate as

[TABLE]

where the unit is in $\mathrm{[kgs^{\bullet-1}]}$ because of the $p$-coordinate.

Let $\tau\in[0,1]$ be a parameter for the order of observational points in a profile. Evaluating the density in Eq. (27) with the state equation $\varrho$, we have

[TABLE]

which has iterated integrals as independent variables. Notice that the shuffle-product property transcribes this as a linear combination of iterated integrals. Substituting this into Eq. (27) finally yields

[TABLE]

This shows that the meridional flow rate $Q_{v}$ is represented as a linear combination of iterated integrals with respect to $T,S,$ and $P$.

Acknowledgements.

The authors appreciate the members of JAMSTEC Argo data management team for preparing and compiling the Argo profile data. All numerical computations were performed on the JAMSTEC DA supercomputer system. Argo float data and metadata are freely available from Global Data Assembly Centre (Coriolis GDAC http://www.coriolis.eu.org/Observing-the-Ocean/ARGO or USGODAE GDAC https://nrlgodae1.nrlmry.navy.mil/argo/argo.html). The processing codes are available on Zenodo [Sugiura (\APACyear2020)].