text:
\begin{document}

\preprint{APS/123-QED}

\title{Sparse Identification of Slow Timescale Dynamics}

\author{Jason J. Bramburger}
\email{jbramburger7@uvic.ca}
\affiliation{Department of Mathematics and Statistics, University of Victoria, Victoria, BC, Canada. V8W 2Y2}
\author{Daniel Dylewsky}
\email{dylewsky@uw.edu}
\affiliation{Department of Physics, University of Washington, Seattle, WA, USA. 98195}
\author{J. Nathan Kutz}
\email{kutz@uw.edu}
\affiliation{Department of Applied Mathematics, University of Washington, Seattle, WA, USA. 98195}

\begin{abstract}
Multiscale phenomena that evolve on multiple distinct timescales are prevalent throughout the sciences.  It is often the case that the governing equations of the persistent and approximately periodic fast scales are prescribed, while the emergent slow scale evolution is unknown. Yet the course-grained, slow scale dynamics is often of greatest interest in practice. In this work we present an accurate and efficient method for extracting the slow timescale dynamics from signals exhibiting multiple timescales that are amenable to averaging. The method relies on tracking the signal at evenly-spaced intervals with length given by the period of the fast timescale, which is discovered using clustering techniques in conjunction with the dynamic mode decomposition. Sparse regression techniques are then used to discover a mapping which describes iterations from one data point to the next. We show that for sufficiently disparate timescales this discovered mapping can be used to discover the continuous-time slow dynamics, thus providing a novel tool for extracting dynamics on multiple timescales.
\end{abstract}

\maketitle

\section{\label{sec:Intro}Introduction}

Many physical phenomena exhibit multiscale dynamics where the fastest timescale is relatively simple to both observe and predict, while the emergent long timescale dynamics are unknown. Examples abound in physics and range from tidal amplitudes \cite{Wunsch,Wunsch2}, to molecular dynamics simulations~\cite{andersen,karplus}, to atmospheric dynamics \cite{Palus}, to the motion of the planets \cite{Brasser,Laskar,Malhotra}. Mathematically, systems which exhibit multiscale dynamics are expensive to simulate since the fastest scales must be accurately resolved, and when the fast timescale is approximately periodic, this expense is primarily used to simulate predictable dynamics. In systems where the scale separation is explicit, it is sometimes the case that one can average the system to obtain the leading-order dynamics of the slow timescale evolution \cite{Guckenheimer,Kifer,Sanders}, although many systems lack an obvious separation of scales or even governing equations, meaning that novel methods for extracting the slow timescale dynamics must be developed.

In this work we present a computationally cheap and efficient method that integrates machine learning and multiscale modeling for extracting and forecasting slow timescale dynamics of a multiscale system. Unlike averaging techniques, the proposed mathematical architecture learns a nonlinear dynamical system characterizing the slow-scale behavior. The setting for the problem is that we are given a signal $\Phiv(t,\varepsilon t)$ in $\mathbb{R}^d$ which contains terms dependent on $t$ contributing to the fast dynamics and terms dependent on $\varepsilon t$, for some $0<\varepsilon \ll 1$, contributing to the slow dynamics. We assume that the fast timescale terms are periodic, so $\Phiv$ satisfies
\begin{equation}\label{Signal}
	\Phiv(t+T,\varepsilon (t+T))=\Phiv(t,\varepsilon (t+T))\approx\Phiv(t,\varepsilon t),
\end{equation}
for all $t$ and some $T > 0$. The timescales of $\Phiv$ are assumed to be sufficiently disparate so that the slow dynamics are nearly constant on the intervals $[t,t+T]$, coming from the size of $\varepsilon$ relative to the fast period $T$. It is through these assumptions that we see that the signal exhibits a fast oscillation and a slow drift, making it amenable to averaging.

Suppose the signal is given on some finite timescale and we wish to understand the physics governing the slow timescale evolution so that the signal can be reconstructed and forecast far beyond the given time window. Since the fast-scale dynamics are relatively simple we would like to `average' these dynamics out to forecast only the slow-scale variable. Knowing the fast-scale period $T>0$ naturally leads to tracking the signal after each period:
\begin{equation}\label{xns}
	\xv_n = \Phiv(0,\varepsilon nT), \quad n = 0,1,2,3,\dots
\end{equation}
since $\Phiv(nT,\cdot) = \Phiv(0,\cdot)$ by periodicity. It follows that understanding the signal at $t = nT$ requires an understanding of the slow timescale dynamics. Our goal is to discover a mapping $\Fv:\mathbb{R}^d \to \mathbb{R}^d$ so that
\begin{equation}\label{SlowMap}
	\xv_{n+1} = \Fv(\xv_n),
\end{equation}
for all $n \geq 0$. Beyond this, if we suppose there exists a function $\Gv:\mathbb{R}^d \to \mathbb{R}^d$ such that the slow timescale dynamics of a signal satisfying Eq.~\eqref{Signal} are governed by the ordinary differential equation (ODE)
\begin{equation}
    \partial_2\Phiv(0,t) = \Gv(\Phiv(0,t)),
\end{equation}
where $\partial_2$ is the derivative with respect to the second component, then expanding the left-hand-side of Eq.~\eqref{SlowMap} as a Taylor series about $t = nT$ using the definition of $\xv_{n+1}$ from Eq.~\eqref{xns} gives that
\begin{equation}\label{Duality}
    \Fv(\xv) = \xv + \varepsilon T \Gv(\xv) + \mathcal{O}(\varepsilon^2).
\end{equation}
Hence, when $\varepsilon$ is small Eq.~\eqref{Duality} represents an Euler step of the slow physics. This demonstrates the duality between the functions $\Fv$ and $\Gv$, at least to $\mathcal{O}(\varepsilon^2)$. Therefore, discovering the mapping $\Fv$ from data presents a novel method to extract the slow timescale physics of a given signal, something which holds great potential for multiscale systems that are so common in the biological and medical sciences \cite{Alber}.

There are potentially two unknowns that are required to find $\Gv$: the fast period $T>0$ and the mapping $\Fv$. Here we use the method of sliding-window {\em dynamic mode decomposition} (DMD)~\cite{Dylewsky} to extract the fast period and the {\em sparse identification of nonlinear dynamics} (SINDy) algorithm \cite{Bramburger,SINDy} to obtain the mapping $\Fv$. We summarize the method visually in Figure~\ref{fig:Slow_Extraction} and explain the individual components in more detail in the following section. The advantage our proposed method presents is that it does not attempt to discover a dynamical system for the multiscale signal, since as we show in Section~\ref{sec:Multiscale}, sparse identification procedures should be expected to fail at this task in general. This is due in part to the numerical approximation of the derivative of the signal, which in the presence of a periodic fast timescale is rapidly changing. Although the challenge of approximating the derivative can be overcome by using instead integral terms \cite{Reinbold,Schaeffer}, the present method circumvents this difficulty altogether. Furthermore, the multiscale property of the data presents a unique advantage in the implementation of a regression technique which promotes sparsity using thresholding: unsupervised timescale separation is expected to introduce error at $\mathcal{O}(\varepsilon^2)$, whereas the sparsification procedure amounts to pruning terms of order $\mathcal{O}(\varepsilon)$. As a result, the method by construction recovers only the leading order slow timescale dynamics.

\begin{figure}[t]

\includegraphics[width=0.5\textwidth]{Slow_Extraction_2.pdf}
\caption{An overview of the method presented in this work. We first use the sliding window DMD technique \cite{Dylewsky} to identify the fast period, then track the signal at integer multiples of the fast period, resulting in the coarsened signal $x_n = \Phiv(0,\varepsilon nT)$. We then identify a mapping $x_n \mapsto x_{n+1}$ using the SINDy method \cite{SINDy}, and finally use the expansion Eq.~\eqref{Duality} to extract the slow timescale physics.}
\label{fig:Slow_Extraction}
\end{figure}

The sparse regression framework has been used previously to discover multiscale physics~\cite{Champion}, but with a much different procedure for the model discovery process and sampling strategy.  Specifically, fast sampling of the dynamics was required since both fast- and slow-timescale physics needed to be discovered, whereas in this work, the fast scale is known and only the slow dynamics needs to be discovered.  Modification to the basic SINDy \cite{SINDy} architecture allows for the discovery of partial differential equations~\cite{Rudy}, equations of rational expressions~\cite{Kaheman}, and conservation laws~\cite{Kaiser}.  For each of these modifications, one can imagine using the proposed method to handle multiscale temporal dynamics in an efficient manner.

This paper is organized as follows. In Section~\ref{sec:Methods} we present a detailed overview of the method for extracting slow timescale physics. This includes summaries of the two major components of the method: the sliding window DMD technique and the SINDy algorithm. In Section~\ref{sec:Multiscale} we provide a simple example of a multiscale system where naive application of the SINDy method fails to properly identify the full governing equations. We further supplement this discussion by showing that our methods accurately capture the slow evolution of this system once the fast dynamics have been taken out of the signal. Section~\ref{sec:Applications} is entirely dedicated to applications of the method, where we present examples with discovered slow dynamics that are monotone, periodic, and chaotic. Finally, we briefly summarize our findings in Section~\ref{sec:Conclusion}.

\section{\label{sec:Methods} Methods}

The method employed by this work has two major components: scale separation using the sliding-window DMD and sparse regression analysis using SINDy. In this section we briefly summarize these methods and direct the reader to the works \cite{Dylewsky} and \cite{SINDy}, respectively, for a more complete discussion of each component. In its totality, our method first employs the sliding window DMD technique to identify the fast period, then track the signal at integer multiples of the fast period, resulting in the coarsened signal $\xv_n = \Phiv(0,\varepsilon nT)$. We then identify a mapping $\Fv$ such that $\xv_{n+1} = \Fv(\xv_n)$ using the SINDy method. Finally, we use the duality Eq.~\eqref{Duality} to extract the slow timescale physics when $\varepsilon$ is sufficiently small. We summarize the method visually in Figure~\ref{fig:Slow_Extraction}, and proceed through the following subsections with a discussion of the individual components of our method.

\subsection{\label{subsec:DMD}Scale separation using sliding-window DMD}

The method presented in this work relies on foreknowledge of the timescale disparity between components constituting a multiscale signal. For systems with known governing equations, this can generally be determined by inspection of the coefficients or by perturbation expansion. For systems with unknown equations of motion, however, a data-driven discovery method must suffice. This could be accomplished by identifying peaks on the Fourier spectrum, for example. For multiscale physics, windowed Fourier transforms can help provide improved resolution of signals and their content.  Indeed, such windowing procedures are the basis of multi-resolution analysis and wavelet decompositions~\cite{mallat}.  Multi-resolution DMD~\cite{mrdmd} provides a analogous decomposition of multivariate data, identifying coherent spatial modes and temporal frequencies in multiscale systems.

Dynamic Mode Decomposition (DMD) is a model regression technique which seeks a best-fit linear representation for observed dynamics~\cite{schmid,dmdbook}. Given a data matrix $\Xv$ consisting of $m$ sequential snapshots, DMD identifies a linear operator $\boldsymbol{A}$ which in some optimal (least-squares) sense satisfies the equation $\dot{\Xv} = \boldsymbol{A}\Xv$. This approximation is of course unlikely to be accurate for highly nonlinear systems, but this problem can be circumvented by subsampling $\Xv$ onto shorter time intervals by sliding a window of width $d$ ($d \ll m$) across the full time series. Even when global dynamics are nonlinear, local linear approximations are easily obtainable. Each local operator $\boldsymbol{A}$ has eigenvalues which characterize the timescale content of local dynamics.

In this work we use a technique introduced by Dylewsky et al. (2019) which leverages DMD to separate the components of multiscale data and learn local linear models for the dynamics by clustering on their eigenspectra \cite{Dylewsky}. This method has the advantages of sparsity, flexibility, and robustness to overlapping or highly disparate dynamical timescales. By gathering eigenspectra across all sliding-window iterations of DMD, one can form a statistical picture of the global spectral content of the data.  A simple clustering algorithm can identify the most prominently represented timescales and offer a parsimonious estimate of the scale components present. Even when timescales are separated by many orders of magnitude, the algorithm can be applied recursively with varied window width to properly identify them. These clusters can be depicted visually by plotting the modulus squared of the eigenvalues of the local operators $\boldsymbol{A}$ for each window. In the present scenario of a signal satisfying Eq.~\eqref{Signal} at least two distinct clusters should be apparent: one bounded away from 0 representing the fast periodic dynamics and another near 0 representing the slow dynamics (which should appear nearly static over the comparatively short span of a windowed subsample). A cartoon of this clustering is presented in panel II of Figure~\ref{fig:Slow_Extraction}, while an application of the method to real data is presented in Figure~\ref{fig:Toy_Cluster}.

\subsection{\label{SINDy}Sparse identification of nonlinear dynamics (SINDy)}

The SINDy method, introduced by Brunton et al. (2016), is an algorithm for discovery of a symbolic representation of the governing equations of a system from time series measurements~\cite{SINDy}. Given a data matrix comprised of sequential state measurement snapshots:
\begin{equation}
    \Xv_1 = \left[
                    \begin{matrix}
                        \text{ ---} & \xv^T(t_1) &  \text{--- }\\
                        \text{ ---} & \xv^T(t_2) &  \text{--- }\\
                         & \vdots &  \\
                        \text{ ---} & \xv^T(t_m) &  \text{--- }
                    \end{matrix}\right]
\end{equation}
along with another data matrix of the same size, $\Xv_2$, comprised of either the temporal derivative of the data at the same measurement times or the successive iterates of the discrete temporal data. A symbolic representation for the dynamics $\dot{\xv}(t) = \boldsymbol{F}(\xv(t))$ or $\xv(t_{n+1}) = \Fv(\xv(t_n))$ is constructed from a library of candidate functions for terms of $\Fv(\xv)$. The chosen functions are evaluated on the measurement data to construct a library matrix $\Theta$ whose $m$ rows represent the $m$ measurement snapshots of $\Xv_1$ lifted into a space of all library observables. For example, a library consisting of polynomials up to degree 2 for $\xv = (x_1,x_2)\in\mathbb{R}^2$ would look like
\begin{equation}
    \begin{split}
    &\Theta(\Xv_1) = \\
    &\begin{bmatrix}
        1 & x_1(t_1) & x_2(t_1) & (x_1(t_1))^2 & x_1(t_1)x_2(t_1) & (x_2(t_1))^2\\
        1 & x_1(t_2) & x_2(t_2) & (x_1(t_2))^2 & x_1(t_2)x_2(t_2) & (x_2(t_2))^2\\
        \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\
        1 & x_1(t_m) & x_2(t_m) & (x_1(t_m))^2 & x_1(t_m)x_2(t_m) & (x_2(t_m))^2\\
      \end{bmatrix}
    \end{split}
\end{equation}
The claim that the equation of motion $\Fv(\xv)$ is some linear combination of the chosen library functions is equivalent to the statement
\begin{equation}\label{eq:sindy_regression}
    \Xv_2 = \Theta(\Xv_1)\Xi \quad
\end{equation}
for some coefficient matrix $\Xi$. For many systems of interest the governing equations contain only a few terms, so a sparseness requirement can be imposed on $\Xi$. This can be carried out using any sparse regression algorithm on the over-determined system of linear equations (\ref{eq:sindy_regression}), e.g. elastic net methods such as LASSO or Ridge. The approach employed in this paper, as in \cite{SINDy}, is a sequentially-thresholded least squares method in which Eq.~\eqref{eq:sindy_regression} is initially solved simply by minimizing $\|\Xv_2 - \Theta(\Xv_1)\Xi\|_2$. Elements of $\Xi$ with magnitude smaller than some chosen sparsity parameter $\lambda > 0$ are then thresholded to $0$, and the process repeats. This iteration is carried out until all remaining coefficients have absolute value larger than $\lambda$. This process is a proxy for $\ell_0$ optimization~\cite{Zheng}. It has convergence guarantees~\cite{Zhang} and performs comparably to LASSO in most cases at significantly reduced computational expense. Moreover, the direct correspondence between the sparsity parameter $\lambda$ and the obtained dynamical coefficients in $\Xi$ proves useful in the multiscale case treated in this paper: we show that knowledge of the timescale components present in a signal can be used to determine {\em a priori} a suitable value for $\lambda$.

We note that there are many variants of sparse regression, all of which typically attempt to approximate a solution to an $NP$-hard, $\ell_0$-norm penalized regression.  Sparsity-promoting methods like the LASSO~\cite{tibshirani1996regression,su2017false} use the $\ell_1$-norm as a proxy for sparsity since tractable computations can be performed.   The iterative least-squares thresholding algorithm of the SINDy algorithm promotes sparsity through a sequential procedure.  Recently, Zhang and Schaeffer~\cite{Zhang} have provided a number of rigorous theoretical results on the behavior and convergence of the SINDy algorithm. Specifically, they prove that the algorithm approximates local minimizers of an unconstrained $\ell_0$-penalized least-squares problem.  This allows them to provide sufficient conditions for general convergence, the rate of convergence, and conditions for one-step recovery. As shown in Champion et al.~\cite{champion2019unified}, the SINDy regression framework does not readily accommodate extensions, additional constraints, or improvements in performance.  Thus the optimization formulation is extended to include additional structure, robustness to outliers, and nonlinear parameter estimation using the {\em sparse relaxed regularized regression} (SR3) approach that uses relaxation and partial minimization~\cite{Zheng}.  Rigorous theoretical bounds are provided for the relaxed formulation which has three advantages over the state-of-the-art sparsity-promoting algorithms: 1) solutions of the relaxed problem are superior with respect to errors, false positives, and conditioning; 2) relaxation allows extremely fast algorithms for both convex and nonconvex formulations; and 3) the methods apply to composite regularizers, essential for total variation (TV) as well as sparsity-promoting formulations using tight frames. Indeed, the SR3 formulation was shown to have superior performance (computational efficiency, higher accuracy, faster convergence rates, and greater flexibility) across a range of regularized regression problems with synthetic and real data, including applications in compressed sensing, LASSO, matrix completion, TV regularization, and group sparsity. Between the three recent papers~\cite{champion2019unified,Zheng,Zhang}, rigorous estimates can be established for this model-discovery paradigm.

\section{\label{sec:Multiscale}SINDy and multiscale systems}

\begin{figure}

\includegraphics[width=0.45\textwidth]{QuasiPeriodic.pdf}
\caption{The first component of the signal $\xv(t)$ in Eq.~\eqref{ToySignal}. There are two distinct timescales and the signal is plotted for slightly more than two periods of the slow timescale. The fast period is approximately $\pi/5$.}
\label{fig:QuasiPeriodic}
\end{figure}

In this section we briefly demonstrate how naively applying the SINDy method to a signal in the form of \eqref{Signal} to discover both the fast and slow timescales simultaneously should be expected to fail, thus necessitating the method described in this manuscript.  Alternatively, one can us fast sampling strategies to resolve the discovery process~\cite{Champion}, but here sampling on the slow scale is all that is needed.  In~\cite{Dylewsky} a simple toy model was created to extract the periods of oscillation for signals formed as linear combinations of periodic phenomena. The model is given by
\begin{equation}\label{ToyModel}
	\begin{split}
		\dot{v}_1 &= v_2, \\
		\dot{v}_2 &= -w_1^2v_1^3, \\
		\dot{w}_1 &= w_2, \\
		\dot{w}_2 &= -100w_1 - 4w_1^3,
	\end{split}
\end{equation}
where constants are chosen to appropriately separate the timescales. The $(w_1,w_2)$ variables are governed by an unforced Duffing equation, for which almost all initial conditions fall into steady periodic motion. The $(v_1,v_2)$ variables form a cubic oscillator with a coefficient $w_1^2$ dependent on the state of $w_1$. The signal \eqref{Signal} is produced by integrating system~\eqref{ToyModel} using MATLAB's ode23 function with a maximal timestep of $10^{-4}$ and initial conditions $(v_1,v_2,w_1,w_2) = (0,0.5,0,0.5)$. We then take a randomly generated orthogonal matrix $Q \in \mathbb{R}^{4\times 4}$ to define
\begin{equation}\label{ToySignal}
	\xv(t) = Q\cdot[v_1(t),v_2(t),w_1(t),w_2(t)]^T,
\end{equation}
in an effort to sufficiently mix the disparate temporal dynamics. We direct the reader to Figure~\ref{fig:QuasiPeriodic} for a visual depiction of the first component of $\xv(t) = (x_1(t),x_2(t),x_3(t),x_4(t))$, and note that the rest look similar in that we can clearly see the two timescales present in the signal $\xv(t)$.

\begin{figure}
\center
\includegraphics[width=0.45\textwidth]{Toy_Comparison.pdf}
\caption{A comparison of the first component of the training data against the evolution of the discovered continuous-time SINDy model with the same initial data. Notice that the discovered data has a significantly higher frequency of oscillation, as can be observed by comparing the extent of the horizontal axis in this figure with that of Figure~\ref{fig:QuasiPeriodic}.}
\label{fig:Toy_Comparison}
\end{figure}

If we use the SINDy method for discovery of a continuous-time dynamical systems using Eq.~\eqref{ToySignal} as the training data and a sparsity parameter $\lambda = 10^{-3}$, we discover the differential equation
\begin{equation}\label{Toy_Discovery}
	\begin{split}
		\dot{x}_1 &= 0.1 - 788.8x_1 + 788.0x_4  + 3.3x_4^3, \\
		\dot{x}_2 &= -833.3x_1 + 833.3x_4, \\
		\dot{x}_3 &= -0.1 - 877.8x_1  + 878.7x_4 - 3.3x_4^3, \\
		\dot{x}_4 &= -0.3 - 922.3x_1  + 924.0x_4 - 6.6x_4^3.
	\end{split}
\end{equation}
Even without knowing the exact values in $Q$ we immediately note that this discovered system cannot exactly capture the dynamics of the signal since it includes constant terms. In Figure~\ref{fig:Toy_Comparison} we plot the first component of the training data against the dynamics of the discovered model \eqref{Toy_Discovery} with the same initial condition. We can see that the discovered model produces a signal with a faster frequency of oscillation than that of the original signal. This can be observed by comparing the extent of the horizontal axes of Figures~\ref{fig:QuasiPeriodic} and \ref{fig:Toy_Comparison}, where the former displays the signal for $t \in [0,100]$ while the latter uses $t \in [0,1]$. This significantly shorter time window in Figure~\ref{fig:Toy_Comparison} is meant to emphasize the extremely fast frequency of oscillation in the discovered continuous-time SINDy model, so it should be noted that the slow dynamics of Eq.~\eqref{ToySignal} are are not yet visible on such a short timescale. We further note that changing the sparsity parameter $\lambda$ has little effect on the discovered model, always producing a dynamical system which fails to even approximately reproduce the dynamics of the input signal \eqref{ToySignal}. The repository {\bf GitHub/jbramburger/Slow-Discovery} contains all code related to this example to ensure the ability to reproduce these results and perform further experiments.

This example is merely meant to illustrate what a naive application of the continuous-time SINDy method can result in when applied to multiscale training data. We suspect that the discrepancy between the SINDy model and the training data is due at least in part to the approximation of the derivative. The fast frequency of oscillation present in the multiscale training data results in a rapidly oscillating derivative, which requires finely spaced data. Hence, one way to improve the performance of the SINDy model might be to take smaller timesteps in the training data, but we note that decreasing the maximal step size down a full order of magnitude to $10^{-5}$ shows little improvement in the discovered continuous-time SINDy model. Of course, in the case of real data it seems unlikely that such a re-generation of the data could be performed and is therefore not a practical method for improving the performance of the discovery algorithm. Another method would be to avoid taking numerical derivatives altogether and instead use weak formulations which replace the derivative with integral terms \cite{Reinbold,Schaeffer}. The strength of the method presented herein is that we are able to consider coarse-grained data which experiences little variation from one data point to the next (representing the slow dynamics), thus avoiding this difficulty altogether.

We now turn to coarsening the signal by tracking it at integer multiples of the fast timescale period. The training data in this case is exactly the same as that used to discovery the continuous-time SINDy model \eqref{Toy_Discovery}. In Figure~\ref{fig:Toy_Cluster} we plot the resulting frequencies $\omega$ extracted from the sliding window DMD method, where one can clearly see the separation between fast and slow dynamics. The centroid of the fast cluster gives that the fast component of the signal has period given by $T = \pi/5$, and hence we can track the full signal at integer multiples of this fast period. We comment that the frequencies showing the largest variation from the clusters are localized to windows near the extreme points of the slow dynamics, representing regions where the slow dynamics are not approximately constant in the relatively small windows. It is exactly this problem that necessitates the sliding window method in the first place, as looking at the whole signal will inevitably introduce large discrepancies since the slow dynamics show large variation on long timescales.

Having now obtained the fast period of oscillation, we are able to apply the sparse identification procedure to obtain a mapping over the slow dynamics. The training data for the SINDy method uses 200 fast periods - long enough to observe at least two full slow periods. The iterates of the component $x_1(nT)$ to the discovered mapping $\Fv$ are presented in Figure~\ref{fig:QuasiPeriodic_Map} against the training data, while we note that the other components look similar. The cumulative error over all components is 1\%, a significant improvement from the model \eqref{Toy_Discovery}. Hence, we see that in the absence of the fast oscillations the SINDy method is able to effectively track the slow evolution of the signal.

\begin{figure}
\center
\includegraphics[width=0.45\textwidth]{om_spec.pdf}
\caption{The result of the sliding window DMD method applied to the signal~\eqref{ToySignal}. The horizontal axis represents the value of $t$ at which the window is centred, while the vertical axis plots the square modulus of the purely imaginary eigenvalues $\omega$ of the local operator {\bf A} for each window. Colours are used to delineate between the two clusters, representing the fast and slow components of the signal. The horizontal black line represents the centroid of the fast frequency cluster, approximately $|\omega|^2 = 100$.}
\label{fig:Toy_Cluster}
\end{figure}

\begin{figure}
\includegraphics[width=0.45\textwidth]{QuasiPeriodic_Map1.pdf}

\caption{Iterates of the discovered mapping for the slow timescale dynamics of the first component of \eqref{ToySignal}, seeded with the initial value of the training data. Training data is provided in blue along with the iterates of the discovered mapping in black.}
\label{fig:QuasiPeriodic_Map}
\end{figure}

\section{\label{sec:Applications}Applications}

Let us consider a simple and motivating example. Consider the scalar ODE
\begin{equation}\label{Logistic}
	\dot{x} = \varepsilon x(1 - x + \sin(2\pi t)), \quad x(0) = x_0 \in\mathbb{R}_+
\end{equation}
with $0 < \varepsilon \ll 1$. Since the ODE is $1$-periodic in the independent variable $t$, the theory of averaging for dynamical systems \cite{Guckenheimer,Sanders} dictates that there exists a constant $C > 0$ such that the solution $x(t)$ satisfies $|x(t) - y(t)| \leq C\varepsilon$ for all $t \geq 0$ \footnote{Technically, the theory dictates that this bound only holds for $t$ on a timescale of length $1/\varepsilon$, but since all positive initial conditions of Eq.~\eqref{Logistic} evolve towards a global attractor, the bound can be shown to hold for all $t \geq 0$ so long as the bound $C > 0$ depends on $x_0$.}, where $y(t)$ is a solution of the autonomous ODE
\begin{equation}\label{AveragedLogistic}
	\dot{y} = \varepsilon y (1-y), \quad y(0) = x_0.
\end{equation}
The solution $y(t)$ is explicitly given by
\[
	y(t) = \frac{x_0}{x_0 + (1 - x_0)\exp(-\varepsilon t)},
\]
where we can therefore see that $y(t)$ evolves on a slow timescale, $\varepsilon t$. Hence, $y(t)$ makes up the slow dynamics to at least $\mathcal{O}(\varepsilon)$ and using Eq.~\eqref{Duality} we expect the mapping \eqref{SlowMap} to be of the form $\Fv(x) = x + \varepsilon x(1-x) + \mathcal{O}(\varepsilon^2)$.

To illustrate the performance of our method, we take $\varepsilon = 10^{-2}$ and a sparsity parameter $\lambda = \varepsilon^2$ to truncate at order $\varepsilon^2$. We discover the mapping \eqref{SlowMap} here to be
\begin{equation}\label{SINDyLogistic}
    \begin{split}
	\Fv(x) &= 1.001x - 0.009778x^2 \\
	&\approx x + \varepsilon x(1-x) + \mathcal{O}(\varepsilon^2),
    \end{split}
\end{equation}
conforming with our expectation from the above analysis. Furthermore, the leading order terms in Eq.~\eqref{SINDyLogistic} represent a forward Euler discretization of Eq.~\eqref{AveragedLogistic} with step size $\varepsilon$, implying that standard error bounding arguments based on the initial condition and the value of $\varepsilon$ can be applied to bound the difference between the solutions of \eqref{Logistic} at integer values of $n\geq 0$ and the iterates of the map \eqref{SINDyLogistic}.

This example illustrates how the method performs against a benchmark ODE where the slow dynamics can be explicitly determined via analysis. In reality, few systems which exhibit multiscale dynamics take the form of singularly perturbed dynamical systems and therefore even determining the fast timescale period to average over is a nontrivial task.

\subsection{\label{subsec:Planets}Planetary dynamics}

As in the toy model example of Section~\ref{sec:Multiscale}, multiscale phenomena can often be observed when there are two periodic components of the signal with a vast separation between their periods. That is, consider a signal which can approximately be written as
\begin{equation}\label{DoublePeriodic}
	\xv(t) = P_0(t) + P_1(t),
\end{equation}
where $P_0$ is periodic with period $T_0> 0$ and $P_1$ is periodic with period $T_1 > 0$. Assuming that $T_0 \ll T_1$, naturally leads to the scale separation parameter
\[
	\varepsilon := \frac{T_0}{T_1} \ll 1.
\]
We may rescale $t = T_0\tau$ so that $\tilde{P}_0(\tau) := P_0(T_0\tau)$ is now $1$-periodic and $P_1(T_0\tau)$ is $\varepsilon^{-1}$-periodic. Setting
\[
	\tilde{P}_1(s) := P_1(\varepsilon^{-1}s)
\]
makes $\tilde{P}_1$ $1$-periodic as well. Hence, the full signal $x(\tau)$ is can equivalently be written
\[
	\tilde{\xv}(\tau) := \xv(T_0\tau) = \tilde{P}_0(\tau) + \tilde{P}_1(\varepsilon \tau)
\]
for which both functions $\tilde{P}_0$ and $\tilde{P}_1$ are $1$-periodic and clearly demonstrate the separation of timescales of interest in this work. The mapping \eqref{SlowMap} tracks $\tilde{\xv}(n)$ for integers $n \geq 0$, or equivalently, $\xv(nT_0)$.

Such a phenomenon was exemplified in the model~\eqref{ToyModel} and can further be found in the motion of the planets Saturn and Jupiter around the sun, as described by a three-body planetary model. For this test case we generate data by constructing a system of three objects subject to pairwise gravitational attraction (using the non-relativistic Newtonian formulation). The masses and initial conditions of these objects are assigned based on measurements of the true Jupiter-Saturn-Sun system in its center-of-mass reference frame. Time series trajectories in $\mathbb{R}^3$ are obtained via symplectic numerical integration over a period of $10^6$ years.

The signal for each planet contains multiple different timescales, and the motion of each planet in its orbital plane can be approximated by a signal of the form~\eqref{DoublePeriodic}. The periodic function $P_0$ describes the primary orbit of the planets around the sun, for which upon applying the sliding window DMD technique we find periods $T_0 = 11.86$ years for Jupiter and $T_0 = 29.5$ years for Saturn. The function $P_1$ describes the eccentricity of these orbits, for which the sliding-window DMD technique has shown that $T_1 \approx 46800$ years~\footnote{It was shown in \cite{Dylewsky} that Jupiter has a periodic component with period between our $T_0$ and $T_1$, but this period is almost exactly $9T_0$ and is expected to be a numerical artifact.} for both planets \cite{Dylewsky}. This periodic eccentricity of the orbits comes from the interaction of the two planets orbiting the relatively massive sun and can be observed in Figure~\ref{fig:Orbits} where we plot their orbits around the sun in their orbital plane.

\begin{figure}
\center
\includegraphics[width=0.5\textwidth]{Orbits.pdf}
\caption{The procession of Jupiter (red) and Saturn (blue) about the sun (black dot) in their orbital plane in nondimensionalized coordinates. The fast timescale for each planet constitutes simple procession about the sun: 11.86 years for Jupiter and 28.5 years for Saturn. Beyond this fast procession are a number of slow scales that drive the eccentricity of the orbit, giving the appearance that they trace out thick circles in the orbital plane over thousands of years.}
\label{fig:Orbits}
\end{figure}

Let us begin by considering our signal to be the motion of Jupiter projected entirely into its orbital plane over time. Here we have $T_0 = 11.86$ with $t$ measured in years and so from the construction above we get $\varepsilon \approx 2.6\times 10^{-4}$. Denote $x^J_n$ and $y^J_n$ to be the $x$ and $y$ components of the position of Jupiter in its orbital plane at time $t = nT_0 = 11.86n$ years. Taking a sparsity parameter $\lambda$ on the order of $\varepsilon^2 \sim 10^{-8}$ is impractical since numerical error present in either simulating the data or implementing the SINDy method should be expected to show up in the discovered mapping. To overcome this, we note that $T_0\varepsilon \sim 10^{-3}$, and so properly discovering $\Gv$ per Eq.~\eqref{Duality} is expected to succeed with $ 10^{-3} \leq \lambda < T_0\varepsilon$ since the lower bound should be large enough to eliminate numerical error and $\mathcal{O}(\varepsilon^2)$ terms, while the upper bound allows for discovery of the $\mathcal{O}(\varepsilon)$ terms. With $\lambda = 10^{-3}$ we discover the slow-scale mapping
\begin{equation}\label{JupiterForecast}
	\begin{bmatrix}
		x^J_{n+1} \\ y^J_{n+1}
	\end{bmatrix} = \begin{bmatrix}
		0.0036 \\ -0.0018
	\end{bmatrix} + \begin{bmatrix}
		0.9999 & 0.0114 \\ -0.0114 & 0.9999
	\end{bmatrix} \cdot \begin{bmatrix}
		x^J_{n} \\ y^J_{n}
	\end{bmatrix}
\end{equation}
using training data comprised of 1500 fast timescale periods - long enough to observe at least two full periods of the first slow timescale. The mapping \eqref{JupiterForecast} can be used to find that the slow timescale mapping $\Gv$ in this case is given by
\begin{equation}\label{JupiterSlow}
    \Gv(x,y) = \begin{bmatrix}
		0 & 1 \\ -1 & 0
	\end{bmatrix}\cdot\begin{bmatrix}
		1.17 + 3.69x \\ 0.57 + 3.69y
	\end{bmatrix}
\end{equation}
using Eq.~\eqref{Duality}. Hence, we see that the eccentricity dynamics of Jupiter is given, to leading order, by an ellipse in the orbital plane. We note that the conservative structure of the slow dynamics is not necessarily guaranteed by the conservative structure of the original three-body problem since it could be the case that energy flows from one scale to the next. Despite this, the resulting conservative structure of Eq.~\eqref{JupiterSlow} could potentially reflect that the training data is stable up to $\mathcal{O}(\varepsilon^2)$ perturbations on very long timescales. Hence, we do not expect that Eq.~\eqref{JupiterSlow} is valid for all time, but only long finite timescales until evolution on even slower timescales begins to influence the $\mathcal{O}(\varepsilon)$ dynamics of Jupiter's orbit.

We may proceed in the same way for the procession of Saturn around the sun to find a similar result. In this case we have $T_0 = 29.5$ years, making up the fast scale period, and therefore the scale separation parameter is given by $\varepsilon \approx 6.3\times 10^{-4} $. Performing our mapping discovery procedure again with $\lambda = 10^{-3}$, for the same reasons as in the case of Jupiter, results in the mapping
\begin{equation}\label{SaturnForecast}
	\begin{bmatrix}
		x^S_{n+1} \\ y^S_{n+1}
	\end{bmatrix} = \begin{bmatrix}
		-0.0038 \\ 0.0025
	\end{bmatrix} + \begin{bmatrix}
		1.0000 & -0.0100 \\ 0.0100 & 1.0001
	\end{bmatrix} \cdot \begin{bmatrix}
		x^S_{n} \\ y^S_{n}
	\end{bmatrix}
\end{equation}
where the variables $(x^S_n,y^S_n)$ represent the orthogonal components of the position of Saturn in its orbital plane at $t = 29.5n$ years. We can similarly rearrange Eq.~\eqref{SaturnForecast} via Eq.~\eqref{Duality} to find that the eccentricity dynamics of Saturn are confined to an ellipse in the orbital plane.

\subsection{\label{subsec:Chaos}Chaotic slow dynamics}

In this example we apply the method to systems for which the slow timescale dynamics are chaotic. We consider a signal
\begin{equation}\label{ChaosSignal}
	\xv(t) = \varepsilon P(t) + C(\varepsilon t),
\end{equation}
where $P(t)$ is periodic with period $T > 0$ and $C(t)$ is a trajectory of a chaotic dynamical system. For the purpose of illustration we will take $C(t)$ to be a trajectory on the circularly symmetric attractor of Thomas~\cite{Thomas}, given by the three-dimensional dynamical system
\begin{equation}\label{Chaos}
	\begin{split}
		\dot{C}_1 &= \sin(C_2) - 0.2C_1, \\
		\dot{C}_2 &= \sin(C_3) - 0.2C_2, \\
		\dot{C}_3 &= \sin(C_1) - 0.2C_3,
	\end{split}
\end{equation}
with initial condition $C(0) = (0.3,0.2,0.1)^T$. The damping value $0.2$ is chosen so that the system is indeed chaotic and the initial condition was chosen arbitrarily to produce a trajectory on the chaotic attractor. The periodic signal $P(t)$ will take the form of a truncated Fourier series
\begin{equation}\label{PerChaos}
	P(t) = a_0 + \sum_{n = 1}^N \bigg[a_n\cos\bigg(\frac{2\pi t}{T}\bigg) + b_n\sin\bigg(\frac{2\pi t}{T}\bigg)\bigg]
\end{equation}
for some $N > 1$ fixed and coefficients $a_0,a_n,b_n\in\mathbb{R}^3$. For our work here we will fix $N = 10$ and note that working with larger or smaller $N$ produces nearly identical results. We further take the coefficients $a_n$ and $b_n$ to be drawn from the uniform distribution on the unit cube $[-1,1]^3$, fix $\varepsilon = 0.1$ and $T = \frac{1}{4}$, which efficiently separates the timescales. The reader is referred to Figure~\ref{fig:Chaos} for a characteristic illustration of the signal.

\begin{figure*}

\includegraphics[width=\textwidth]{Thomas1.pdf} \\
\includegraphics[width=\textwidth]{Thomas2.pdf}
\caption{A signal of the form \eqref{ChaosSignal} with periodic fast dynamics given by Eq.~\eqref{PerChaos} with randomly-generated Fourier coefficients and slow dynamics governed by the system~\eqref{Chaos}. The top shows the signal over a long timescale so that the chaotic dynamics can be observed. The bottom shows a smaller time window so that only the periodic fast scale can be observed.}
\label{fig:Chaos}
\end{figure*}

As before we will fix $\lambda = \varepsilon^2 = 0.01$ to ensure that only $\mathcal{O}(\varepsilon)$ terms are present in the discovered equation. However, since the discovered equation contains sinusoidal functions, a library containing only monomials will never be able to fully reproduce the dynamics of the slow scale. Hence, the best one could hope for is to produce a monomial series representation of the dynamics of \eqref{Chaos}, but the truncation at $\mathcal{O}(\varepsilon^2)$ will inevitably truncate the series representation as well. Furthermore, the chaotic nature of system~\eqref{Chaos} leads one to conjecture that SINDy will not necessarily just return a truncated Taylor series representation for $\sin(x)$ since the elements of $C(t)$ do not remain small on large timescales. This is confirmed by implementing the discovery method 100 different times, resulting in 100 different randomized functions $P(t)$, and discovering $\Fv(\xv) = \xv + \mathcal{O}(\varepsilon^2)$ in every case, meaning that the slow timescale dynamics are unable to be picked up at all with such a monomial basis. This emphasizes that the choice of basis functions plays potentially an even more critical role in the case of discovering the dynamics of the slow timescale system than it potentially would for discovering systems with a single timescale.

The inclusion of both sine and cosine terms into the library immediately remedies the above problem. For example, one implementation of this method resulted in the system
\begin{equation}\label{DiscoverChaos}
	\Fv(\xv) = \begin{bmatrix}
		0.02508\sin(x_2) + 0.9949x_1  \\
		0.02428\sin(x_3) + 0.9947x_2  \\
		0.02471\sin(x_1) + 0.9951x_3
	\end{bmatrix},
\end{equation}
which is characteristic of all implementations with randomized Fourier coefficients in Eq.~\eqref{PerChaos}. After rearranging, the slow dynamical system $G(\xv)$ is recovered as
\[
    \Gv(\xv) = \begin{bmatrix}
	    1.0032\sin(x_2) - 0.2040x_1  \\
		0.9712\sin(x_3) - 0.2120x_2 \\
		0.9884\sin(x_1) - 0.1960x_3
	\end{bmatrix},
\]
which agrees with Eq.~\eqref{Chaos} up to $\mathcal{O}(\varepsilon^2)$. \vspace{.1 in}

\section{\label{sec:Conclusion}Conclusion}

In this work we have seen a computationally cheap and efficient method for discovering the slow timescale dynamical system of a multiscale signal that is amenable to averaging. This method relies on tracking the signal not continuously, but after each fast period, to produce a mapping which parsimoniously describes the slow physics. The resulting mapping can also be related to the continuous-time dynamics of the slow timescale via the Euler stepping scheme Eq.~\eqref{Duality}. Furthermore, the fast period can be extracted accurately using the recently developed sliding window DMD technique with clustering of eigenfrequencies \cite{Dylewsky}, which can also be used to reconstruct the fast component of the signal. The result is an algorithm that reliably reproduces the slow evolution of a multiscale signal using DMD for the fast timescale and SINDy for the slow timescale. Given the difficulty in approximating the emergent slow scale evolution dynamics of multiscale systems, the method provides a viable architecture for coarse-graining to achieve accurate, interpretable and parsimonious dynamical models for slow-scale physics. Finally, we note that depending on the scaling, averaging can yield trivial slow dynamics, in which case homogenization is the correct procedure \cite{Givon}. This potentially limits the class of signals for which our techniques are applicable, although our work herein provides an initial step towards data-driven discovery of the slow timescale dynamics of multiscale signals.

\begin{acknowledgments}
JJB was supported by a PIMS PDF held at the University of Victoria. JNK acknowledges support from the Air Force Office of Scientific Research (AFOSR) grant FA9550-17-1-0329.
\end{acknowledgments}


TITLE: Magic squares using matrices
QUESTION [0 upvotes]: Let $M(d)_{d\in \mathbb{R}}\subset M_3(\mathbb{R})$, consisting of all the matrices such that the sum of the elements in each row, column and on the two diagonals is $d$. Such an element $M(d)$ is called a magic square with magic number $d$.
Show that $M(0)$ is a vector subspace of $M_3(\mathbb{R})$ and find a basis for $M(0)$. Conclude that the center entry of any element in $M(0)$ is $0$.
Is there another way to doing this problem other than guessing the form of $M(0)$? I thought about writing down $\left(\begin{array}{ccc}
\lambda_1 & \lambda_2 & \lambda_3 \\
\lambda_4 & \lambda_5 & \lambda_6 \\
\lambda_7 & \lambda_8 & \lambda_9 \end{array}\right)$ and summing up the rows and columns to zero, but seems unproductive

REPLY [1 votes]: The sum of the first and third rows plus the sum of the first and third columns, minus the sum of the two diagonals,  gives
$$
2\lambda_2 + 2\lambda_8 + 2\lambda_4 + 2 \lambda_6 -2 \lambda_5 = 0.
$$
Similarly, the sum of the middle row and middle column is zero, but it's also
$$
4 \lambda_5 + Q
$$
where $Q$ is the sum of the four "mid-edge" entries (2, 4, 6, 8).
Taking twice the second equation minus the first gives
$$
8 \lambda_5 + 2Q - (2Q - 2 \lambda_5) = 0 \\
10 \lambda_5 = 0.
$$


TITLE: Example of context-free grammar that triggers exponential behaviour without memoization in RD parsers
QUESTION [4 upvotes]: It is often said that memoization brings the complexity of recursive-descent parsers
from exponential to polynomial. However, I had a hard time finding an example grammar
that triggers the exponential case.
Because recursive descent parsers cannot (naturally) deal with left-recursion,
the example should not contain left-recursive rules.
BTW, I'm familiar with PEG and Packrat parsers, and I'm not looking for examples
in the that domain. To be more precise, can someone please give me a non-left-recursive
context-free grammar that triggers exponential behaviour?

REPLY [7 votes]: Here is a grammar that should meet your specification, though it
generates the very simple language $a^+(b+c)$.  (A simpler grammar has been added below)
$S \rightarrow ab \mid aBb \mid ac \mid aCc$
$B \rightarrow a \mid aB \mid aBB$
$C \rightarrow a \mid aC \mid aCC$
How was it built:
I started from the grammar
$\{S \rightarrow S S \mid a\}$ which is very ambiguous, so that trying all
parses of a string $a^n$ takes exponential time.
I put it in Greibach Normal Form, so as to remove any left recursion,
giving
$S \rightarrow a \mid aB$
$B \rightarrow a \mid aB \mid aBB$
Then I introduced $C$ to duplicate the possibilities of derivation,
with $b$ or $c$ as end marker to discover at the very end whether I
should have been using $B$ or $C$ as non-terminal in my derivation,
i.e,, whether the first rule to apply was $S \rightarrow aBb$ or $S
\rightarrow aCc$. But if the choice of rule made at the very beginning
was wrong, the backtracking process will have first to try all parsing
possibilities for the string of $a$, which has exponential cost.
Since the recursive descent has to try first one of the two rules,
there will always be a string (requiring the other rule) that takes
exponential time to parse, because of the useless backtracking with the wrong non-terminal.
Memoisation brings this down to cubic time. Naive memoisation will make it $O(n^4)$, but a bit of currification can correct that to cubic time.

A simpler solution (thought of later)
All that matters in the above grammar is to force an early choice of
the terminal marker $b$ or $c$, while having that choice checked only at the end
of the parsing process.
Hence the following grammar is enough to get the desired exponential
behaviour. Actually it is a (renaming) homomorphic image of the first
grammar above, and leads to exactly the same computation up to that
renaming.
$S \rightarrow ab \mid aXb \mid ac \mid aXc$
$X \rightarrow a \mid aX \mid aXX$
The rules $S \rightarrow ab \mid ac$ are not actually necessary for
our purpose, but that is minor. Removing them would change the language to $aa^+(b+c)$. Then, it is also possible to remove the initial $a$ in the rules $S \rightarrow aXb \mid aXc$ without introducing left recursion, thus getting back the initial language $a^+(b+c)$. But the grammar is no longer in Greibach normal form (which was not a requirement of the problem). This gives a solution with only 5 rules:
$S \rightarrow Xb \mid Xc$
$X \rightarrow a \mid aX \mid aXX$


TITLE: Word probability problem
QUESTION [0 upvotes]: I am trying to solve this probability problem:
Suppose there are 100 students total in your year at UCLA. Each tutorial group consists of 10 students.
What is the probability that for the next course you are in a tutorial group with none of the students
that were in your tutorial group for Basic Statistics (assuming, of course, that students are assigned
to groups in a random manner)?
Now when I read it I think that it could be solved by just calculating the probability complement of odds of meeting 10 people out of 100 which would be:
1 -10/100 = 0.9
However, this question is from a material which discusses binomial distribution, and combinations and permutations, so I feel like I miss something. Any ideas?

REPLY [1 votes]: There are ${99\choose 9}$ ways to assign studuents (other than your self) to join you in any study session.
And ${90\choose 9}$ that include none who are in your stats study session.
$\frac {{90\choose 9}}{{99\choose 9}}$


TITLE: If $T:V\longrightarrow V$ is diagonalizable, then $T$ has a characteristic polynomial of the form $f(x)=(x-c_{1})^{n_{1}}\dotsc(x-c_{k})^{n_{k}}$.
QUESTION [1 upvotes]: If $c_{1},\dotsc, c_{k}$ are the distinct eigenvalues of $T$ and $W_{i}$ is the eigenspace of $c_{i}$, then $T$ has characteristic polynomial of the form $f(x)=(x-c_{1})^{\dim W_{1}}\dotsc (x-c_{k})^{\dim W_{k}}$.
How do I prove this?

REPLY [1 votes]: I found an answer. If $B$ is an eigenbasis of $T$ for $V$, then $[T]_{B}$ is diagonal. Hence for each eigenvalue $c_{i}$, $[T-c_{i}I]_{B}$ is diagonal, and the nullity of a diagonal matrix is equal to the number of zeroes on the diagonal.


TITLE: Taking suitable odd power of $\phi$, to get order of $\phi$ a power of $2$?
QUESTION [0 upvotes]: Let $\phi\in Aut(G)$, for a finite group $G$. We also have $\phi^2\in Inn(G)$. Therefore it is clear that if we can prove that $\phi^m\in Inn(G)$ for any odd $m$, then because $\phi^2$ is inner, it implies that $\phi$ is also inner.
Now in a  proof I was doing, author wanted to assume order of $\phi$ to be a power of $2$ for some purpose. Thus he says by taking suitable odd power of $\phi$, we may assume that order of $\phi$ is a power of $2$.
Now suppose I have $\phi$ which has order 10. Now I need to replace my $\phi$ with another $\phi'=\phi^k$ for some odd $k$, such that order of $\phi'$ becomes a power of $2$. Suppose I take $\phi'=\phi^5$, then I am done.
But now if I have $\phi$ which has order $7$, then what do I take my $\phi'$. Why does this work for all orders of $\phi$.

REPLY [1 votes]: If order of $\phi$ is already odd like $7$, notice that $<\phi>=<\phi^2>$. Hence
$\phi^2\in Inn(G)\implies \phi\in Inn(G)$. Thus, you can assume that order of $\phi$ is an even number.


TITLE: On quantities with no very small odd prime factors; a response to Wlodzimierz Holsztynski
QUESTION [2 upvotes]: In response to a comment posted under
 Powers of $2$ and the products of initial odd primes , I shall raise some questions about quantities near $O_n= P_{n+1}/2$, the product of the first $n$ odd primes. (Note the change in meaning of $O_n$; $P_n$ is the $n$th primorial.)
Before going into the questions, note that for integers $d$ of absolute value less than the largest prime factor of
$O_n$ and coprime to $O_n$, such integers d have only 2 as a prime factor and $O_n + d$ has all its prime factors greater than the $n+1$st prime, unless $d$ has absolute value 1, in which case 2 is also a prime factor.  Let this set of integers be called $C_n$.

(From the other question) How large a power of 2 divides $O_n^2 - 1$?  As seen elsewhere, this quantity is not a power of 2 when $n>1$.  More interestingly, the greater the power of 2 which divides this quantity, the further away one can demonstrate easily $O_n$ is from a power of 2.
Let $f(n)$ be that integer $e$ such that  $2^e$ precisely divides
$O_n^2 - 1$.  So $f(1)=3, f(2)=5, f(3)=4.$ Is $f$ unbounded as a function
of $n$?  How fast does it grow?  Is $f$ the sum of two nonconstant integer valued functions, one of them periodic and the other monotonic?
How many primes are in $C_n$?
Form $Q_n$, the product of the integers in $C_n$.
This product is a multiple of 8, so is not squarefree.  How far is it from being squarefree?  In particular, is $Q_n$ the multiple of the square of any odd prime?
Given recent work on prime gaps (Zhang et al.), can we hope for a Generalized Dirichlet theorem? Specifically, translate the constellation $C_n$ by multiples of $O_n$: are there infinitely many translates which contain two or more primes in that translate?

I hope this makes my comment regarding related questions appear less vague and less broad.

REPLY [2 votes]: Answer to question 2: obviously $f(n)\geq 3$ for all $n$, but I will show that $f(n)=3$ for infinitely many $n$, so we cannot give any other lower bound for the growth of $f$. To see this, notice that if $p_2...p_{n+1}\not \equiv \pm 1 \pmod 8$ (where $p_i$ is the $i$th prime), then the highest power of $2$ in $O_n^2-1$ is $2^3$. On the other hand, if $p_2...p_{n+1}\equiv\pm 1 \pmod 8$ and $n$ is chosen so that $p_{n+2}\equiv \pm 3 \pmod 8$ (by Dirichlet's theorem, there are infinitely many such $n$), we get $p_2...p_{n+2}\not \equiv \pm 1 \pmod 8$, so in any case there are infinitely many $n$ for which $f(n)=3$. We can even say something about the density of such $n$.  By denoting $A=\{n:f(n)>3\}$, $B=\{n:f(n+1)=3\}$ and $C=\{n:p_{n+2}\equiv \pm 3 \pmod 8\}$, the natural denisties satisfy  $d(A)+d(B)=1$ and we have $A\cap C \subset B$. Hence, using the PNT in aritmetic progressions, which tells that $d(C)=\frac{1}{2}$, we get $d(B)\geq d(A)-\frac{1}{2}=\frac{1}{2}-d(B)$, and from this we solve $d(B)\geq \frac{1}{4}$.
Answer to question 3: To clarify the definition,  $C_n=\{O_n\pm d: d<p_{n+2} \hspace{0.1 cm}\text{is a power of two}\}.$ So the number of elements up to $x$ in $\bigcup_{n\geq 1} C_n$ grows like $\log^2 x$. As far as I know, no "natural" set that contains logarithmically few elements has been proved to contain infinitely many primes, so estimating the number of primes in $C_n$ seems hopeless. But one would expect that even $O_n+2$ is prime infinitely often.
Answer to question 4: first notice that the greatest common divisor of any two elements of $C_n$ is one. Indeed, if $q>2$ is a prime satisfying $q\mid O_n\pm 2^{\alpha}, q\mid O_n\pm 2^{\beta}$ and $2^{\alpha}\leq 2^{\beta}<p_{n+2}$, then $q\mid 2^{\beta-\alpha}-1,$ but this is a contradiction since $q\geq p_{n+2}$. So the question of $Q_n$ being divisible by an odd prime square reduces to an element of $C_n$ having this property. Again, this seems extremely hard since for example the following problem seems easier but is known to be open: Is $2^{p}-1$ ($p$ prime) always squarefree?
Answer to question 5: This would at least follow from Dickson's conjecture. According to it, there should even exist translates of $C_n$ that consist entirely  (except $O_n\pm 1$) of primes. The constellation is admissible for this conjecture since if $q>2$ was a prime that divided one element of  $kO_n+C_n$ for every $k$, then $q$ would be coprime to $O_n$ and bigger than $p_{n+1}$, but then $kO_n \pmod q$ would attain all the values $1,...,q-1$ as $k$ varies, which is a contradiction.


TITLE: Show that the inequality $\frac{s}{t} + \frac{t}{r} + \frac{r}{s} >= 3$, assuming s,t,r >0.
QUESTION [0 upvotes]: For this question, I made an assumption that $r<=s<=t$. If I manage to proof it under this assumption, is the result still hold true when $r<=t<=s$ or $ t<=s<=r$, if it is true, why?
Sorry, I didn't learn about AM-GM inequality. Would you mind to explain how AM-GM inequality works in this question here? My knowledge of AM-GM is limited to $a^2+b^2>=2ab $....
Thank you very much for your reply. This question has been confusing me for a long period of time.

REPLY [2 votes]: $$\sum_{cyc}\left(\frac{s}{t}-\frac{s}{r}\right)=\frac{\sum\limits_{cyc}(s^2r-s^2t)}{str}=\frac{(s-t)(t-r)(r-s)}{str}\geq0$$ for $r\leq s\leq t$, which says that it's enough to prove our inequality for $r\leq t\leq s$.
But it's better to use AM-GM: $$\sum_{cyc}\frac{s}{t}\geq3\sqrt[3]{\frac{s}{t}\cdot\frac{t}{r}\cdot\frac{r}{s}}=3.$$
I use AM-GM for three numbers:

for positives $x$, $y$ and $z$ we have:
  $$\frac{x+y+z}{3}\geq\sqrt[3]{xyz}$$

or $$x+y+z\geq3\sqrt[3]{xyz}.$$
Now, take $x=\frac{s}{t},$ $y=\frac{t}{r}$ and $z=\frac{r}{s}.$


text:
\begin{document}

\maketitle

\begin{abstract}

We investigate stochastic optimization problems under relaxed assumptions on the distribution of noise that are motivated by empirical observations in neural network training. Standard results on optimal convergence rates for stochastic optimization assume either there exists a uniform bound on the moments of the gradient noise, or that the noise decays as the algorithm progresses. These assumptions do not match the empirical behavior of optimization algorithms used in neural network training where the noise level in stochastic gradients could even increase with time. We address this behavior by studying convergence rates of stochastic gradient methods subject to changing second moment (or variance) of the stochastic oracle as the iterations progress. When the variation in the noise is known, we show that it is always beneficial to adapt the step-size and exploit the noise variability. When the noise statistics are unknown, we obtain similar improvements by developing an online estimator of the noise level, thereby recovering close variants of RMSProp~\cite{tieleman2012lecture}. Consequently, our results reveal an important scenario where adaptive stepsize methods outperform SGD.

\end{abstract}

\vspace{-0.3cm}
\input{sections/introduction}

\input{sections/motivation}

\input{sections/baseline}
\vspace{-0.2cm}

\input{sections/online_estimation}
\vspace{-0.3cm}

\input{sections/experiment}
\vspace{-0.3cm}

\section{Conclusions}
\vspace{-0.2cm}
This paper discusses convergence rates of stochastic gradient methods in an empirically motivated setting where the noise level is changing over iterations. We show that under mild assumptions, one can achieve faster convergence than the fixed step SGD by a factor that is  polynomial in number of iterations, by applying online noise estimation and adaptive step sizes. Our analysis, therefore provides one explanation for the recent success of adaptive methods in neural network training.

There is much more to be done along the line of non-stationary stochastic optimization. Under our current analysis, there is a gap between the adaptive method and the idealized method when the noise variation is large (see second row in Table~\ref{tab:complexity}). A natural question to ask is whether one could reduce this gap, or alternatively, is there any threshold preventing the adaptive method from getting arbitrarily close to the idealized baseline? Moreover, could one attain further acceleration by combining momentum or coordinate-wise update techniques? Answering these questions would provide more insight and lead to a better understanding of adaptive methods.

Perhaps a more fundamental question is regarding the iterate dependency: the setting where the moments~$m_k$ or the variance~$\sigma_k$ are functions of the current update $x_k$, not just of the iteration index $k$. Significant effort needs to be spent to address this additional correlation under appropriate regularity conditions. We believe our work lays the foundation to address this challenging research problem.

\section{Acknowledgement}

This work was partially supported by the MIT-IBM Watson AI Lab. SS and JZ also acknowledge support from NSF CAREER grant Number 1846088.

{

  \bibliographystyle{abbrvnat}
  \bibliography{bibliography}
}

\newpage

\appendix

\input{sections/appendix}

\end{document}

config: arxiv


TITLE: If $x=y^2+z^2$, $y=z^2+x^2$, and $z=x^2+y^2$, prove $\frac{x}{x+1} + \frac{y}{y+1} + \frac{z}{z+1} = 1$
QUESTION [1 upvotes]: If $x=y^2+z^2$, $y=z^2+x^2$, and $z=x^2+y^2$,
prove
$$\frac{x}{x+1} + \frac{y}{y+1} + \frac{z}{z+1}  = 1$$

REPLY [0 votes]: Hint:
From the first two,
$$x-y=y^2-x^2$$ factors as $$(x-y)(x+y+1)=0$$
which defines a pair of planes, and similarly
$$(z-x)(z+x+1)=0.$$
Eliminate $y,z$ and plug in the last identity.

Note that
$$\frac{-(x+1)}{-(x+1)+1}=\frac{x+1}{x}$$
and the last identity has one of the forms
$$3t=1,\\2t+t^{-1}=1,\\t+2t^{-1}=1.$$


TITLE: Reducing download time using prime numbers
QUESTION [38 upvotes]: An issue has cropped up recently in programming with which I could greatly benefit from the expertise of proper mathematicians. The real-world problem is that apps often need to download huge chunks of data from a server, like videos and images, and users might face the issue of not having great connectivity (say 3G) or they might be on an expensive data plan.
Instead of downloading a whole file though, I've been trying to prove that it's possible to instead just download a kind of 'reflection' of it and then using the powerful computing of the smartphone accurately reconstruct the file locally using probability.
The way it works is like this:

A file is broken into its bits (1,0,0,1 etc) and laid out in a predetermined pattern in the shape of a cube. Like going from front-to-back, left-to-right and then down a row, until complete. The pattern doesn't matter, as long as it can be reversed afterwards.
To reduce file size, instead of requesting the whole cube of data (the equivalent of downloading the whole file), we only download 3 x 2D sides instead. These sides I'm calling the reflection for want of a better term. They have the contents of the cube mapped onto them.
The reflection is then used to reconstruct the 3D cube of data using probability - kind of like asking a computer to do a huge three-dimensional Sudoku. Creating the reflection and reconstructing the data from it are computationally heavy, and as much as computers love doing math, I'd like to lighten their load a bit.

The way I'm picturing is like a 10x10 transparent Rubik's cube. On three of the sides, light is shone through each row. Each cell it travels through has a predetermined value and is either on or off (either binary 1 or 0). If it's on, it magnifies the light by its value. If it is off, it does nothing. All we see is the final amount of light exiting each row, after travelling through the cube from the other side. Using the three sides, the computer needs to determine what values (either 1 or 0) are inside the cube.
At the moment I'm using normal prime numbers as the cell values to reduce computing time, but I'm curious to know if there is another type of prime (or another type of number completely) that might be more effective. I'm looking for a series of values that has the lowest possible combination of components from within that series.
Here is a rough diagram:

It might help to imagine that light shines in at the green arrows, and exits with some value at the red arrows. This happens for each row, in each of the three directions. We're left with only three 2D sides of numbers, which are used to reconstruct what's inside the cube.
If you look where the 14 exits on the left, it can have two possible combinations, (3 + 11 and 2 + 5 + 7). If for arguments sake we were to assume it were 3 and 11 (coloured green), then we could say at the coordinate where 3 and 11 exist, there are active cells (magnifying the light by their value). In terms of data, we would say this is on (binary 1).
In practice we can rarely say for certain (for 2 and 3 we could) what an inside value has based on its reflection on the surface, so a probability for each is assigned to that coordinate or cell. Some numbers will never be reflected on the surface, like 1, 4 or 6, since they can't be composed of only primes.
The same happens in the vertical direction, where the output is 30, which has multiple possibilities of which two correspond to the possibility shown in the horizontal direction with an exit of 14, coloured blue and pink (since they hit the 23, the same as 3 in the horizontal direction).
This probability is also added to that coordinate and we repeat in the front-to-back direction, doing the same a final time. After this is done for each cell in the whole cube, we have a set of three probabilities that a cell is either on or off. That is then used as the starting point to see if we can complete the cube. If it doesn't 'unlock' we try a different combination of probabilities and so forth until we have solved the 3D Sudoku.

The final step of the method is once the cube is solved, the binary information is pulled out of it and arranged in the reverse pattern to how it was laid out on the server. This can then be cut up (say for multiple images) or used to create a video or something. You could cough in theory cough download something like Avatar 3D (280GB) in around 3 minutes on decent wifi. Solving it (nevermind building the pixel output) would take a while though, and this is where I'm curious about using an alternative to prime numbers.
You might have guessed that my maths ability goes off a very steep cliff beyond routine programming stuff. There are three areas of concern / drawbacks to this method:

it is rubbish at low levels of data transfer. A 10 x 10 x 10 cube for instance has a larger 'surface area' than volume. That's because while each cell can hold one bit (either 1 or 0), each surface cell needs to be a minimum of 8 bits (one character is one byte, or 8 bits). We can't even have 'nothing', since we need null to behave as a type of placeholder to keep the structure intact. This also accounts for why in the above diagram, a 1000x1000x1000-cell cube has its surface areas multiplied by 4 characters (the thousandth prime is 7919 - 4 characters) and the 10'000(cubed)-cell cube has its surface areas multiplied by 6 characters (10'000th prime is 104729, six characters). The aim is to keep total character length on the 2D side to a minimum. Using letters could work, as we could go from a-Z with 52 symbols, before paying double bubble for the next character (the equivalent to "10"). There are 256 unique ASCII characters, so that's the upper limit there.
the factorials are still too high using prime numbers. Is there a series of numbers that are both short in character length (to avoid the problem above) and have very few possible parents? I'm leaning towards some subset of primes, but lack the maths to know which - some sort of inverted Fibonacci? The fewest possible combinations, the faster the computer will solve the cube.
I haven't tested yet if its possible to use a third, fourth or nth side to increase either the capacity of the cube or the accuracy of the reflection. Using a say octahedron (yellow below) instead of a cube might be better, just hurts the brain a little to picture how it might work. I'm guessing it would tend towards a sphere, but that's beyond my ability.

EDIT:
Thank you for your helpful input. Many answers have referred to the Pigeonhole principle and the issue of too many possible combinations. I should also correct that a 1000 x cube would require 7 digits not 4 as I stated above, since the sum of the first 1000 primes is 3682913. I should also emphasise the idea isn't to compress in the common sense of the word, as in taking pixels out of an image, but more like sending blueprints on how to build something, and relying only on the computation and knowledge at the receiving end to fill in the blanks. There is more than one correct answer, and will mark correct the one with the most votes. Many thanks for the detailed and patient explanations.

REPLY [1 votes]: It's a comment, just does not fit
As one of those guys with slow connections, appreciate the efforts.
Youtube is good for me exactly because it has 144p videos, and if something is interesting or I need it I can choose which xxxp to set.
But it has its problems switching and buffering which it probably would not have if that variety of resolutions would be implemented/achieved differently.
In the good old days of modems, it was more often seen, a picture which starts as few huge pixels and then detail'es as long as more data is downloaded. Not sure which compression it was, but it looks like a jpeg on steroids. Those data chunks could be defined as n-streams for video or multiple images, and it could give a finer control on the user side, but setting a limit on how many channels to request, or adjust it automatically.
So maybe take look not in the direction of how much to compress, but how to control quality. All fixed bitrate and variable bitrate compression formats do that varying quality or attempting to keep quality. But they produce one data blob, and manage different resolutions may be challenging on the server-side, so the user does not have that many options (glad there are some places where one has), lags switching quality, etc, etc.
So maybe, as a suggestion, approach the problem from a different perspective -more flexible control of quality on user and server-side. More seamless transition from one quality to another and alike.
People on slower connections already do compromises and sure they appreciate places that offer suitable options. Same way as there is no alternative for youtube for me because the least quality I have seen was 360p and it's top-quality for me (with some difficulties, if weather is right so to speak) if they(alternatives) have such options at all. Sometimes I would not mind 100p or do not need video at all until I do, and when I do, I would like it not to discard the already existing 1min buffer but add to it to improve quality.
So better tools and options and approaches to fight all those compromises may be more appreciated than an amazing algorithm that does 10% better than all existing algorithms on any data. By varying quality I get a 90% reduction of data most of the time - when details begin to be important then yes it's time for heavy guns like waiting for a buffer and other available options.
In general, trading user resources for less data transfer is not impossible, and it is not necessarily purely a math problem - like 5 sentences describing an apple tree or a picture, and user side ai generates one or finds visually similar from a dictionary, etc. But it is a bit more vast problem.


TITLE: $M = \{\alpha \in C(I_{b},\overline{B_{a}(0)}):\alpha(0)=x \}$ is closed
QUESTION [1 upvotes]: The following problem has popped up in my study of differential equations.
Consider $C(I_{b},\overline{B_{a}(0)})$, the set of all continuous maps from the compact interval $I_{b} = [-b,b]$ to $\overline{B_{a}(0)}$, the closed ball of radius $a$ centered around $0$. Such maps (i.e. elements of $C$) are obviously bounded, since they are continuous functions on a compact interval.
But consider the subset of $C(I_{b},\overline{B_{a}(0)})$ defined by
$M = \{\alpha \in C(I_{b},\overline{B_{a}(0)}):\alpha(0)=x \}$.
Is this set closed with respect to the norm defined by
$d(f,g) = \sup_{t \in I_{b}}\|f(t)-g(t)\|$
If so, why?

REPLY [1 votes]: Let me write $\mathcal{C} =  C([-b,b],\overline{B}(0,a))$. The evaluation at $0$ function, $\text{ev}_0: \mathcal{C} → \Bbb R$, defined by $\text{ev}_0(\alpha):= \alpha(0)$ is continuous wrt $d$. The set you are looking for is
$$ M = \text{ev}_0{}^{-1}(\{x\}) $$
which is a pre-image of a closed set under a continuous function.


TITLE: Principal ideal $(a) = R$ iff $a$ is a unit
QUESTION [2 upvotes]: I was wondering if for a ring $R$ with $1$, the following is true:
The ideal $(a)$ generated by some element $a \in R$ is the whole ring ($(a) = R$) iff $a \in R^\times$.
The direction $a$ is unit $\Longrightarrow (a) = R$ is clear.
But what about the other direction? I suspect that it might also hold but can't quite prove it.
My attempt is to show that if a non-unit generates $R$, then I need to be able to combine a unit (or in particular $1$) just from $a$ which (hopefully) might not be possibe. Here I get stuck.
Can someone tell me if the implication $(a) = R \implies a \in R^\times$ is correct and help me to prove it, or give me a counterexmaple?
Thanks!
EDIT: To summarize the comments: $R$ has to be commutative otherwise my statement is not true. If $R$ commutative, then the direction I was asking about holds trivially as $(a) = aR$ and $1 \in (a)$ thus $1 = ar$ and $a \in R^\times$.

REPLY [2 votes]: If $(a)=R$, then $1\in (a)$. In particular, $ra=1$ for some $r\in R$, and thus $a$ is a unit.


TITLE: Minimal polynomial over an extension field divides the minimal polynomial over the base field
QUESTION [9 upvotes]: I need help proving this theorem:

Given the field extension: $\mathbf{K} \subseteq \mathbf{L}$, for $\alpha \in \mathbf{L}$ and $g(x) \in \mathbf{K}[x]$, $\alpha$'s minimal polynomial over $K$,
  and $f(x) \in \mathbf{L}[x]$, $\alpha$'s minimal polynomial over $L$,
  then the degree of $g$ is bigger than the degree of $f$ and $f(x)$ divides $g(x)$.

REPLY [5 votes]: Because $\mathbf{K}\subseteq\mathbf{L}$, you also have $\mathbf{K}[x]\subseteq\mathbf{L}[x]$, so that $g\in \mathbf{L}[x]$ and $g(\alpha)=0$, and therefore (because $f$ is the minimal polynomial of $\alpha$ over $\mathbf{L}$) we must have $f\mid g$, and hence also $\deg(f)\leq\deg(g)$.


TITLE: Prove that $\lim_{(x,y) \to (0,0)} \frac{\sin^2 xy}{x^2 + y^2}=0$ without using inequalities
QUESTION [2 upvotes]: I read a similar question that was solved by using the fact that $|\sin(x)| \leq |x|$, but I tried a different approach, as I struggle with utilizing inequalities to solve limits of this kind. Please note that my textbook did not ask for the limit to be solved without using inequalities.
I attempted to prove this result by first switching to polar coordinates, then applying Hopital's rule as follows:
$$\lim_{(x,y) \to (0,0)} \frac{\sin^2 xy}{x^2 + y^2}=\lim_{\rho \to 0} \frac{\sin^2 (\rho^2\sin(\theta)\cos(\theta))}{\rho^2}$$
Now I use cosine duplication formula to rewrite $\sin^2(t)$
$$\lim_{\rho \to 0} \frac{\sin^2 (\rho^2\sin(\theta)\cos(\theta))}{\rho^2}=\lim_{\rho \to 0} \biggl(\frac{1}{2}\biggr)\frac{1-\cos(2\rho^2\sin(\theta)\cos(\theta)}{\rho^2}$$
and now I apply Hopital's rule
$$\lim_{\rho \to 0} \biggl(\frac{1}{2}\biggr)\frac{4\rho\sin(\theta)\cos(\theta)\sin(2\rho^2\sin(\theta)\cos(\theta))}{2\rho}=\lim_{\rho \to 0} \sin(\theta)\cos(\theta)\sin(2\rho^2\sin(\theta)\cos(\theta))$$
Now if I'm understanding this correctly, because the argument of the second sine function goes to zero, the result is proven.
Have I made any mistakes? Was there a faster or more intuitive approach to solving this problem without using inequalities?

REPLY [3 votes]: Another way
First rewrite the expression:
$$
\frac{\sin^2 xy}{x^2 + y^2}
= \left( \frac{\sin xy}{xy} \right)^2 \frac{x^2y^2}{x^2+y^2}
$$
Here the first factor tends to $1$ since $\frac{\sin t}{t} \to 1$ as $t \to 0$ and $xy \to 0$ as $(x,y) \to (0,0).$
For the second factor we use polar coordinates:
$$
\frac{x^2y^2}{x^2+y^2}
= \frac{\rho^4 \cos^2\theta \sin^2\theta}{\rho^2}
= \rho^2 \cos^2\theta \sin^2\theta
\to 0
$$
as $\rho \to 0$ which is the case when $(x,y) \to (0,0).$
Thus,
$$
\frac{\sin^2 xy}{x^2 + y^2}
\to 1^2 \cdot 0
= 0.
$$


TITLE: Finding the conjugates, why can they argue this way?(exercise)
QUESTION [1 upvotes]: In one exercise I am supposed to find the conjugate of $\sqrt{2}+i$ over $\mathbb{Q}$.
I found the answer by finding irr$( \sqrt{2}+i,\mathbb{Q})$, and then solving the polynomial finding all the roots(this is the definition of conjugates), this gave the correct answer, but in the solutions manual, they use a trick I do not understand.
They say that:
The conjugate are $\sqrt{2}+i,\sqrt{2}-i,-\sqrt{2}+i,-\sqrt{2}-i$. This is clear beacuse $\mathbb{Q}(\sqrt{2}+i)=(\mathbb{Q}(\sqrt{2}))(i)$.
Can someone say how this is clear and why it follows from the fact that: $\mathbb{Q}(\sqrt{2}+i)=(\mathbb{Q}(\sqrt{2}))(i)$?

REPLY [1 votes]: By what you mention in your comments, the conjugates of $(\sqrt{2}+i)$ can be computed if we knew the image of $(\sqrt{2}+i)$ under every homomorphism
$$
\varphi : \mathbb{Q}(\sqrt{2}+i) \to \overline{\mathbb{Q}}
$$
Since
$$
\mathbb{Q}(\sqrt{2}+i) = \mathbb{Q}(\sqrt{2})(i)
$$
we need to determine $\varphi(i)$ and $\varphi(\sqrt{2})$. However, since $\varphi$ is a field homomorphism
$$
\varphi(\sqrt{2})^2 = \varphi(2) = 2
$$
since any field homomorphism must fix $\mathbb{Q}$ (this takes a short proof). However, $2$ has exactly 2 square roots in $\overline{\mathbb{Q}} \subset \mathbb{C}$, so
$$
\varphi(\sqrt{2}) = \pm \sqrt{2}
$$
Similarly,
$$
\varphi(i) = \pm i
$$
and so the list of possible conjugates is
$$
\{\pm\sqrt{2}\pm i\}
$$
To prove that each of these is, in fact, a conjugate of $(\sqrt{2}+i)$, you need the fact that the number of such homomorphisms is precisely equal to the degree
$$
[\mathbb{Q}(\sqrt{2}+i):\mathbb{Q}]
$$
My answer is long enough already, so let me know if this helps and if you have any questions.


TITLE: How to work out family of functions graphs?
QUESTION [0 upvotes]: I would like to know manually how to draw the family functions.
I try to find the pictures of the graph, apparently i can find only $x^3$ graph photo below.
what I want to know is how to work out and get this graph?
i have questions like  function is $f(x)=x^n$
1. $x^1$
2. $x^{-1}$
3. $x^{-1/4}$
4. $x^{1/4}$
5. $x^{-2/3}$
but I would like to know how to draw the graph by manually solve the functions.
Hope my question is clear.
Thanks in advance.
eg.  $x^3$

REPLY [1 votes]: You can plot functions using Wolfram|Alpha, e.g. $x^{-2/3}$. If you want to do it by hand, make a table of some function values in the range of interest, mark the corresponding points in a coordinate system and try to connect them with a smooth curve, using whatever you may now about the function, e.g. zeros, poles, extrema etc.


TITLE: How do I determine equilibration in a Monte Carlo $NVT $ simulation?
QUESTION [15 upvotes]: I'm running an NVT (constant number of particles, volume and temperature) Monte Carlo simulation (Metropolis algorithm) of particles in two dimensions interacting via Lennard-Jonse potential $$U = 4\left(\frac{1}{r^{12}} - \frac{1}{r^6}\right),$$ (in reduced units). Boundary conditions are periodic.
From this simulation I'm calculating the instantaneous pressure and potential energy. In the first steps the system is not in equilibrium, so I need to start averaging after the system is in equilibrium.
I'm starting my simulation from a random configuration.
My question: even after the system has reached equilibrium, it fluctuates around this equilibrium. These fluctuations may be large for large temperatures. So how do I know that I have reached equilibrium?
Here are some examples of the curve:

$$\uparrow$$ The energy Vs. simulation step, for a high temperature (warmer color is higher density)

$$\uparrow$$The energy Vs. simulation step, for a low temperature (warmer color is higher density)

$$\uparrow$$The energy Vs. simulation step, for a high temperature, only for low densities. in this graph it's harder to tell if we reached equilibrium (warmer color is higher density)

REPLY [0 votes]: You can't know for sure.
In the simulation, the system is constantly changing chaotically and in general it is possible that after some period of relative constancy of macroscopic variables, some radical change such as phase transition will occur.
However, if the system is very simple, such as hard sphere gas, one can calculate theoretically how the equilibrium state should look like for given constraints such as volume and temperature. Then, if the values derived from the simulation approach those theoretical values, we can say that most probably we are observing a transition to equilibrium and nothing surprising is going to happen.
If we do not know the values of variables in thermodynamic equilibrium, we can't say the system is approaching it.
However, there is also the idea of restricted thermodynamic equilibrium, defined only for some limited time period that we are interested in. The system is then in equilibrium, if estimation of average kinetic energy is uniform throughout the system's region of space and if its observed fluctuation in time and space is consistent with Boltzmann's distribution, and if the same is true for other physical quantities of interest (average pressure, net potential energy).


TITLE: Internal energy and difference in volume at constant temperature
QUESTION [1 upvotes]: If the volume of a cylinder of a gas is decreased so that the work done on the gas is $W$ and the temperature is kept constant by taking away heat energy $H$.
$\Delta U = -H + W$.
Someone told me that if the temperature of the gas is constant does $\Delta U=0$ Which means $W=H$.
But I don't agree to this because

$U$ is sum of all kinetic AND potential energies of the gas. When volume of the gas decrease the average separation of the molecules decrease which means that average potential energy of the molecules will decrease.
As temperature, $T \propto E_{k}  $ so $E_{k} $ is constant.
$E_{k} $ is constant but potential energy decrease which means $U$ change.

So is my conjecture wrong?
Note: i am talking about real gases with significant attractive forces.

REPLY [2 votes]: Who is right depends on what your model for the gas is.  In short, you are right for models that incorporate inter-particle interactions.
The internal energy of an ideal gas, which is a model that neglects interactions between particles, can be written only in terms of the number of particles $N$ in the gas sample and the absolute temperature $T$ of the gas:
\begin{align}
  E = \alpha NkT,
\end{align}
where $\alpha$ is a dimensionless constant.  So in the case of an ideal gas, a sample with a constant particle number and temperature will indeed have $\Delta U = 0$.  Since in introductory courses one usually only considers ideal gases, confusion often results because particularly astute students such as yourself bring up the precise objections you have.
If one models the gas by including inter-particle interactions, then the above argument fails.  Consider, for example, the Van der Waals gas.  Its internal energy depends not only on $N$ and $T$, but also on $V$, the volume of the gas, precisely for the reason you outlined.


TITLE: Exponential and log functions compose to identity
QUESTION [5 upvotes]: How to prove that the exponential function and the logarithm function are the inverses of each other? I want it the following way. We must use the definition as power series, and must verify that all the terms of the composition except the coefficent of $z$ vanish, and that the first degree term is $1$.
I can write down the proof for the coefficient of $x^n$ for arbitrary but fixed $n$ by explicit verification. But how to settle this for all $n$ at one go?

REPLY [2 votes]: Qiaochu Yuan's proof is without doubt the most elegant. However, here is a purely formal argument using the derivative:
Let $L(x) = \log \frac 1{1-x} = \sum_{n\geq 1} \frac{x^n}n$. I will show that $\exp(L(x))=\frac 1{1-x}$.
Let $u(x)=\exp(L(x))$. By the chain rule for differentiation, we have
$$
u'(x) = L'(x)u(x).
$$
Noting that $L'(x)=(1-x)^{-1}$, the above can be rewritten as
$$
u(x)=(1-x)u'(x).
$$
Suppose $u(x)=\sum_{n\geq 0} a_nx^n$, then comparing the coefficients of each power of $x$ in the preceding identity gives
$$
a_n = (n+1)a_{n+1} - na_n
$$ for each $n\geq 0$.
Also, $a_0=1$.
It follows that $a_n=1$ for all $n$, so $u(x)=\frac 1{1-x}$, as desired.
The reverse identity, $\log(\exp(x))=x$ is simpler; the chain rule tells us that the derivative of $\log(\exp(x))=1$.


TITLE: Alternative to Ziegler's "Lectures on Polytopes"
QUESTION [8 upvotes]: I am interested in alternatives to Ziegler's Lectures on Polytopes, which is the suggested textbook for a course I am attending. I find the conversational style of the book jarring.

REPLY [2 votes]: Some of the topics of Ziegler's (excellent) book (e.g. Gale transforms, f-vectors, the secondary polytope, fiber polytopes) are also covered in the "Triangulations" book of De Loera, Rambau and Santos. The treatment is from a slightly different perspective and definitely worth a look. They provide many pictures.
Another alternative (also more from the triangulations perspective) is the book Lectures in Geometric Combinatorics.


TITLE: show that $\bigcap_{n=1}^{\infty}{A_n}=A_1 \setminus \left( \bigcup_{n=1}^{\infty} \left( A_1 \setminus A_n\right)\right)$.
QUESTION [0 upvotes]: I don't have advances, I can't understand... I need a hint please!
If $A_n \subseteq X$ for $n=1,2,3,\ldots$,
show that
$$\bigcap_{n=1}^{\infty}{A_n}=A_1 \setminus  \left( \bigcup_{n=1}^{\infty}
 \left( A_1 \setminus A_n \right)\right)\;.$$

REPLY [1 votes]: For those who (like me) prefer a more calculational proof style, and equivalence to mutual implication, here is an alternative 'element chasing' proof, as a very late answer.$%
\require{begingroup}
\begingroup
\newcommand{\calc}{\begin{align} \quad &}
\newcommand{\op}[1]{\\ #1 \quad & \quad \unicode{x201c}}
\newcommand{\hints}[1]{\mbox{#1} \\ \quad & \quad \phantom{\unicode{x201c}} }
\newcommand{\hint}[1]{\mbox{#1} \unicode{x201d} \\ \quad & }
\newcommand{\endcalc}{\end{align}}
\newcommand{\Ref}[1]{\text{(#1)}}
%$
Letting $\;n\;$ range over the positive integers, we calculate which $\;x\;$ are in the right hand side set:
$$\calc
    x \in A_1 \setminus \bigcup_n (A_1 \setminus A_n)
\op\equiv\hint{definitions of $\;\setminus\;$, $\;\bigcup\;$, and $\;\setminus\;$}
    x \in A_1 \;\land\; \lnot \langle \exists n :: x \in A_1 \land x \not\in A_n \rangle
\op\equiv\hint{logic: assume LHS of leftmost $\;\land\;$ in its RHS; simplify}
    x \in A_1 \;\land\; \lnot \langle \exists n :: x \not\in A_n \rangle
\op\equiv\hint{logic: DeMorgan}
    x \in A_1 \;\land\; \langle \forall n :: x \in A_n \rangle
\op\equiv\hint{RHS implies LHS}
    \langle \forall n :: x \in A_n \rangle
\op\equiv\hint{definition of $\;\bigcap\;$}
    x \in \bigcap_n A_n
\endcalc$$
By set extensionality, this proves the original statement.
$%
\endgroup
%$


TITLE: Simple Inverse Matrix Inequality
QUESTION [2 upvotes]: I believe this inequality is true, but I am not sure how to show it rigorously.
For Matrices $A^{-1}$ and $B^{-1}$ where $A^{-1}, B^{-1} \preceq I$ where I is a matrix of all 1's and this inequality (I believe) represents that each element is less than 1, I want to show the following
$$
||A^{-1} - B^{-1}|| = ||A^{-1}A^{-1}A - B^{-1}B^{-1}B|| \le ||A-B||
$$
My reasoning behind me doing this is that I know, $|1/5-1/3| \le |5-3|$ and we can replace 5 or 3 with any natural numbers and this will work.
So how can I rigorously show this for my matrix inequality above? Or is it actually fine as is, as long as I state $A^{-1}, B^{-1} \preceq I$ ?
Thank you for taking the time!

REPLY [0 votes]: Use identity $A^{-1}-B^{-1}=A^{-1}(B-A)B^{-1}$ and the fact that $\|A^{-1}\|\leq 1$ and $\|B^{-1}\|<1$.


TITLE: Evaluating a definite integral
QUESTION [1 upvotes]: Evaluate the definite integral. $$\int_3^7 \ln(x^{43})\,dx$$

REPLY [2 votes]: This is sort of a joke answer (but correct).
If we do not know any "log laws," let $u=\log(x^{43})$, $dv=dx$. Then
$du=\dfrac{43x^{42}}{x^{43}}\,dx=\dfrac{43}{x}\,dx$, and we can take $v=x$.
Thus our indefinite integral is
$$x\log(x^{43})-\int \frac{43}{x} x\,dx.$$


TITLE: 4 points not separable by SVM
QUESTION [0 upvotes]: We know in a support vector machine:

Considering we have a linear feature mapping $\phi(x_n)=x_n$ and the XOR problem.
We have 2 classes in $R^2$, class 1 $ t_+=+1$  and class 2 $t_-=-1$ and 4 points where $x_1, x_2$ are from class1 and $x_3, x_4$ are from class2.
Therefore we can write, $w=a_1x_1+a_2x_2-a_3x_3-a_4x_4$ and $a_1+a_2-a_3-a_4=0$
How can we prove the 4 points are not separable?
Coordinates for $x_1=1,1; x_2=-1,-1; x_3=1,-1, x_4=-1,-1$?

REPLY [0 votes]: the straight forward thing to do it is to assume that it is possible. Note that separability is something that is independent from the SVM. so we just have to proof, that there is no separating hyperplane and don't need to find an SVM solution. Thus we can optimize in the space of $w$.
This gives the following 4 inequalities the solution must satisfy:
$w^T x_1 +b > 0\\
w^T x_2 +b > 0\\
w^T x_3 +b < 0\\
w^T x_4 +b < 0$
now insert the coordinates of the points and try to solve for $w_1$ and $w_2$. You will see, that this is not possible.