text:
\section{#1}\vspace{0pt}}
\newcommand{\mysubsection}[1]{\vspace{0mm}\subsection{#1}\vspace{0mm}}
\DeclareMathOperator{\spn}{span}

\def\sectionautorefname~#1\null{%
  Section~#1\null
}
\def\subsectionautorefname~#1\null{%
  Subsection~#1\null
}
\def\equationautorefname~#1\null{%
  (#1)\null
}

\newcommand{\myfig}[2][placeholder]	%
{
  \begin{figure}[htb]
    \centering
    \ifthenelse{\boolean{singlecolumn}}
    {
      \begin{minipage}[c]{.5\textwidth}
	\includegraphics[width=\textwidth]{figs/#1}
      \end{minipage}\hfill
    }
    {\includegraphics[width=\linewidth]{figs/#1}}
    \ifthenelse{\boolean{singlecolumn}}
    {
      \begin{minipage}[c]{.4\textwidth}
      }
      {\vspace{-8mm}}
      \caption{#2}
    \label{fig:#1}
      \ifthenelse{\boolean{singlecolumn}}
      {
      \end{minipage}
    }
    {}
  \end{figure}
}

\newcommand{\myfigfull}[2][placeholder]	%
{
  \begin{figure*}[htb]
    \centering
    \includegraphics[width=\textwidth]{figs/#1}
    \caption{#2}
    \label{fig:#1}
  \end{figure*}
}

\newenvironment{nalign}{
    \begin{equation}
    \begin{aligned}
}{
    \end{aligned}
    \end{equation}
    \ignorespacesafterend
}

\bibliographystyle{ieeetr}

\section{MRI Physics}\label{sec:physics}

\subsection{Spin Dynamics}\revnum{3}{1}
The MRI system can be approximated by a dynamical system based on the Bloch equations \cite{bib:bloch1946qs}, introduced by Felix Bloch and illustrated conceptually in \autoref{fig:mri_dynamical_system}. An aggregate of spins in each voxel creates a net magnetization that is initially aligned with the scanner's main magnetic field, and evolves based on intrinsic biophysical tissue parameters and user control inputs consisting of radio-frequency (RF) pulses and magnetic field gradients. RF pulses act to rotate the magnetization vector away from the main (longitudinal) field direction and toward the transverse plane. The strength and duration of the RF pulse determine the degree of rotation, or \textit{flip angle}, experienced. The transverse component of the magnetization is sensed by nearby receive coils through Faraday's Law of Induction. We will denote the transverse magnetization as $M_{xy}$ and the longitudinal magnetization as $M_z$.
\myfig[mri_dynamical_system]{The MRI dynamical system, governed by the Bloch equations. An aggregate of spins at each voxel position creates a net magnetization, initially pointed in the longitudinal direction, that evolves based on user pulse sequence control inputs and intrinsic tissue parameters. The acquired image is the transverse component of the magnetization.}

Although many tissue parameters influence the signal evolution, here we limit our scope to relaxation parameters, which are the most common contrast mechanism used in MRI.
Relaxation is a fundamental component of nuclear magnetic resonance and dictates the rate that magnetization returns to the equilibrium state, effectively ``resetting'' the MRI system. Longitudinal magnetization exponentially recovers to its initial state with time constant $T_1$\xspace, and transverse magnetization exponentially decays to zero with time constant $T_2$\xspace.

As \autoref{fig:Bloch_Toilet} shows, relaxation can be understood through an analogy with a toilet\footnote{Introduced by Al Macovski in the 2009 ISMRM Lauterbur Lecture.}, where the water in the tank represents longitudinal magnetization, water in the bowl represents transverse magnetization, and the toilet flush is an RF excitation.
When the toilet is flushed (excitation), water transfers from the tank to the bowl, producing a detectable signal. As the toilet bowl drains ($T_2$\xspace relaxation), the tank refills with water ($T_1$\xspace recovery). Successive flushes will transfer new water from the tank to the bowl.
Neglecting other system effects, the magnetization evolution after a 90$^\circ$ RF pulse is given by
\begin{align}
  M_z(t) &= 1 - e^{-\frac{t}{T_1}}, \label{eqn:t1} \\
  M_{xy}(t) &=  e^{-\frac{t}{T_2}} \label{eqn:t2},
\end{align}
where $M_{z}(t)$ and $M_{xy}(t)$ are the longitudinal and transverse magnetization components at time $t$ after the excitation, respectively, and have initial conditions of $M_z(0) = 0$ and $M_{xy}(0) = 1$ immediately following the RF pulse.\revnum{3}{2}

\myfig[Bloch_Toilet]{Signal relaxation visualized through the toilet analogy. The water in the tank represents longitudinal magnetization, water in the bowl represents transverse magnetization, and the toilet flush is an RF excitation. When the toilet is flushed (excitation), water transfers from the tank to the bowl, producing a detectable signal. As the toilet bowl drains ($T_2$\xspace relaxation), the tank refills with water ($T_1$\xspace recovery). Successive flushes will transfer new water from the tank to the bowl.}

In general, the transverse magnetization distribution in space and time can be described by
\begin{align}
  M_{xy}(\vec r, t) = \rho(\vec r)f_t\left(\vec \theta(\vec r), \vec u_\tau(\vec r)|_{\tau=0}^t\right),
  \label{eqn:mri:mxyt}
\end{align}
where $\vec r \in\real{3}$ represents the 3D spatial position, $\rho(\vec r)$ is the intrinsic amount of magnetization at position $\vec r$ (called proton-density, or PD), and $f_t(\cdot)$ is a spatio-temporal signal evolution that depends on both the biophysical tissue parameters given by $\vec \theta(\vec r)$ and the full history of user-controllable pulse sequence and scanner parameters\footnote{Although the pulse sequence parameters may not have a spatial component, the scanner hardware may introduce a spatial dependence, e.g. due to spatial inhomogeneities.} given by $\vec u_\tau(\vec r)$, $0 \leq \tau \leq t$.\revnum{1}{3}
Although there are numerous tissue parameters and user controls that influence the magnetization, we will focus on the case of FSE imaging, in which the magnetization is sampled at $T$ equispaced intervals with spacing $T_s$. In this case, the image of interest at position $\vec r$ and state $i$ is the sampled transverse magnetization and is given by
\begin{align}
    x_i(\vec r) = M_{xy}(\vec r, i T_s) = \rho(\vec r)f_{iT_s}\left(\vec \theta(\vec r), \vec u_\tau(\vec r)|_{\tau=0}^{iT_s}\right).
    \label{eqn:fse_signal_evolution}
\end{align}

For FSE, the magnetization is primarily sensitive to relaxation parameters, i.e., $\vec \theta = \left( T_1, T_2 \right)$, and to RF refocusing flip angles, i.e. $\vec u_{iT_s}(\vec r) = \mathrm{RF}_i(\vec r)$, where $\mathrm{RF}_i$ represents the flip angle of the $\nth{i}$ RF pulse.\revnum{1}{4}
Although the RF pulse is not prescribed for each position, in practice the flip angles smoothly vary spatially due to RF transmit field inhomogeneity effects that impart varying levels of RF power across the imaging volume.
We represent the vector of magnetization points at the echo times by $\vec f\left(\vec \theta(\vec r), \vec u(\vec r)\right) \in \complex{T}$, where the $\nth{i}$ component of $\vec f$ is equal to $f_{iT_s}$ (i.e. the transverse magnetization signal).\revnum{1}{5}
Based on the sequence timing and RF flip angle inputs, different types of image contrasts can be created. \autoref{fig:brain_multicontrast} shows four common FSE image contrasts primarily due to PD, $T_1$\xspace, and $T_2$\xspace, and created by using different sequence parameters in independent scans.
\myfig[brain_multicontrast]{Different image contrasts based on PD, $T_1$\xspace, and $T_2$\xspace produced by careful choice of RF flip angles and sequence timing for FSE-based scans.}

In particular, $\vec f\left(\vec \theta(\vec r), \vec u(\vec r)\right)$ can be modeled by solutions to the Bloch equations
\cite{bib:bloch1946qs}, which are differential equations that describe the magnetization evolution as a function of time\footnote{In fact the Bloch equations themselves are also phenomenological!}. The Bloch equations allow us to calculate the magnetization (signal) evolution for individual spins in a spatial region, given pulse sequence inputs. For example, the relaxation behavior given by \autoref{eqn:t1} and \autoref{eqn:t2} is the solution to the Bloch equations when no time-varying fields are present. Many extensions have been introduced to model additional contrast mechanisms, including diffusion %
and chemical exchange, %
though these are beyond our scope.

\subsection{Simulating Spins}
Simulating the magnetization evolutions of spins can be a powerful tool for developing and evaluating physics constrained reconstruction methods. For example, in Section~\ref{sec:relax_mri}, we will show how to use simulated signal evolutions to construct an approximate linear subspace for reconstruction. The Bloch equations are linear differential equations, and their numerical simulation can be efficiently computed under reasonable assumptions and approximations.\revnum{3}{4}
In particular, when the RF pulses are discretized, spin simulation using the Bloch equations reduces to successive applications of rotation followed by relaxation to the magnetization vector $\vec{M} = (M_{xy}, M_z)$. Figure~\ref{fig:simulation} shows an illustration of the simulation process. The overall signal evolution is non-linear but differentiable with respect to the relaxation and system parameters.\revnum{2}{3}

\myfigfull[simulation]{When the RF pulses are discretized, spin simulation using the Bloch equations reduces to successive applications of rotation followed by relaxation to the magnetization vector $\vec{M}$. An alternative to the Bloch equations is the extended phase graph (EPG) formalism, which simulates the signal evolution of a distribution of spins across a voxel. EPG efficiently keeps track of signal evolutions across multiple resonant frequencies using the Fourier series.}

In real imaging applications, multiple resonant frequencies often appear within each voxel due to local magnetic field inhomogeneities.\revnum{1}{6} That is, the magnetic field experienced by spins within a voxel can be slightly different. For a realistic simulation, the Bloch equations have to be solved for each resonant frequency separately that appears within one particular voxel, and then summed together to form the resulting signal. This can be computationally demanding.

An alternative to the Bloch equations is the extended phase graph (EPG) formalism, which simulates the signal evolution on a voxel level~\cite{bib:weigel2015}. EPG uniformly discretizes the resonant frequencies within one voxel. The discretization allows EPG to efficiently keep track of signal evolutions across multiple resonant frequencies using the Fourier series. Similar to the Bloch equations, spin simulation using the EPG involves successive applications of rotation followed by relaxation to the underlying magnetization representation. The overall signal evolution using EPG is also non-linear but differentiable.

\section{MRI Sampling}
\subsection{Spatial Encoding}
The received MRI signal represents spatial frequencies of the transverse magnetization distribution in space, and the linear relationship is described by the integral
\begin{align}
  s(t) = \int_{\vec r} M_{xy}(\vec r, t) e^{-j 2 \pi \vec{k}(t)^\top\vec r} d\vec{r} + w(t),
  \label{eqn:signal_equation}
\end{align}
where $s(t) \in \mathbb{C}$ is the acquired signal at time $t$, $\vec{k}(t) \in\real{3}$ is a trajectory through the 3D frequency space, and $w(t)$ is complex-valued white Gaussian noise. The symbol $(\cdot)^\top$ denotes the transpose operation and $j=\sqrt{-1}.$
Because the spatial frequency wave-number is typically denoted as $\vec k$, MRI acquisitions are often described as sampling in \textit{k-space} \cite{bib:fessler_model_based}.

Throughout this article we consider a discrete Fourier approximation of \autoref{eqn:signal_equation}.
Though clinical MRI systems come standard with parallel imaging receive arrays \cite{bib:yingspm}, we will omit their discussion for brevity, with the understanding that they can be flexibly incorporated into the sampling model.
We first consider the acquisition model when the magnetization stays constant over time, i.e. $M_{xy}(\vec r, t) = M_{xy}(\vec r)$ with no relaxation effects. Then, given an underlying image $\vec x \in \complex{N}$, which consists of the transverse magnetization of all $N$ voxels, the full forward model is represented in matrix form as
\begin{align}
    \vec y = \vec P \vec F \vec x  + \vec w,
    \label{eqn:forward_model_multicoil_full}
\end{align}
where $\vec y\in\complex{M}$ are the acquired k-space measurements, $\vec F\in\complex{N\times N}$ is a discrete Fourier transform operator, $\vec P \in \complex{M\times N}$ is a sampling operator that selects the acquired k-space measurements, and $\vec w \in \complex{M}$ is the noise.
The encoding operator is succintly represented as $\vec E = \vec P \vec F$.

Since the scan time is directly proportional to the number of measurements, we are typically interested in solving problems for the case where $M < N$. Compressed sensing offers an avenue for targeting this regime by exploiting low-dimensional structure in the image representation \cite{bib:lustig2008pe}. A common inverse problem approach to CS MRI is a regularized least-squares optimization given by
\begin{align}
\underset{\vec x}{\arg\min}\frac{1}{2}\norm{\vec y - \vec E \vec x}_2^2 + \lambda \mathcal{R}(\vec x),
  \label{eqn:lasso}
\end{align}
where $\mathcal{R}$ is a sparsity-promoting regularization, e.g.\ $\ell_1$ norm of  the wavelet coefficients, and $\lambda>0$ is a regularization term.

\subsection{Multi-dimensional Extensions}
Although initial work on CS MRI focused on sparsity of static, anatomical images using spatial wavelet transforms and total variation \cite{bib:lustig2008pe}, many extensions have been proposed to handle additional imaging dimensions, including joint sparsity, low rank, and their variants \cite{bib:trzasko2013dw, bib:frankmlr, bib:christodoulou2018}.
In particular, the linear forward model can be extended to represent additional image states, $\vec x = \begin{bmatrix}\vec x_1 & \cdots & \vec x_T\end{bmatrix}$, where $\vec x_i \in \complex{N}$ is the image at the $\nth{i}$ image state and $T$ is the number of states. The forward model, illustrated in \autoref{fig:forward_model_SA} for signal relaxation, includes different encoding operators for each state $i$
based on the user-specified sampling patterns:
\begin{align}
  \vec E_{i} = \vec P_i \vec F.
  \label{eqn:temporal_forward_model}
\end{align}
This  concept can be used to include additional dimensions representing signal relaxation \cite{bib:t2shmethod}, cardiac and respiratory motion \cite{bib:lustig2008pe,bib:feng2016zk,bib:jiang2018pe}, and many others \cite{bib:christodoulou2018}. Selecting the $T$ sampling patterns in a way that maintains compatibility with parallel imaging and CS is an active research area \cite{bib:levine2018jk, bib:haldar2018nj}.
  \myfig[forward_model_SA]{The forward model extended to temporal relaxation. Each image represents a sampling time $i$ during the acquisition with $T$ time points, where each time point is Fourier transformed and sampled with a different sampling operator $\vec P_i$ (represented by red circles).}

\section{Model-based MRI}\label{sec:model_based_mri}
The signal evolution model developed in the previous sections describes how the signal received by the MRI scanner is formed. Combining with spatial encoding, the physics based forward model considers the following non-linear evolution:
\begin{nalign}
    x_i(\vec r) &= \rho(\vec r)f_{iT_s}\left(\vec \theta(\vec r), \vec u(\vec r)\right),\\
    \vec{y}_i &= \vec{E}_i \vec{x}_i, \quad i=1,\dots,T.
    \label{eqn:physics_based_forward_model}
\end{nalign}
The most explicit use of the physics model in a reconstruction is to directly solve for the tissue parameters from the raw k-space measurements. This is in contrast to first reconstructing a time series of images followed by a parameter fit. It is
possible to formulate this model-based inversion even in the case of undersampled k-space  \cite{bib:block2009,bib:zhao2014b}.
This can be written as a non-linear, non-convex least squares objective in which we aim to solve for ($\vec \rho$, $\bm \theta$):
\begin{nalign}
    & \underset{\vec \rho, \bm \theta}{\text{minimize}}
& & \frac{1}{2} \norm{\vec{E} \vec{x} - \vec{y}}_2^2 + \lambda \mathcal{R}(\vec \rho, \bm \theta)\\
& \text{subject to}
& & x_i(\vec r) = \rho(\vec r)f_{iT_s}\left(\vec \theta(\vec r), \vec u(\vec r)\right),\\
& & & i=1,\dots,T.
\label{eqn:model_based_recon}
\end{nalign}
In the MRI literature, equation~\autoref{eqn:model_based_recon} is often referred to as quantitative MRI with a model-based reconstruction, as the physical model is incorporated into the objective function.

Compared to the conventional image series reconstruction followed by a fit, solving for the parameters directly serves as the ultimate dimensionality reduction: instead of computing a set of images, the aim is to recover only the intrinsic information represented by $\vec \rho$ and $\vec \theta$. The main benefit of this approach is that the problem size is significantly reduced. For example, for mapping PD and $T_2$\xspace, solving a problem of the form \autoref{eqn:lasso} consists of $TN$ unknown variables, while the model-based representation \autoref{eqn:model_based_recon} has only $2N$ unknowns. In addition, sparsity-promoting regularization penalties can be applied directly to the parameter maps.
Since the data fidelity term is differentiable, the overall problem can be optimized using first-order or second-order methods.

Additional system-related parameters such as RF field inhomogeneity and off-resonance can also be modeled and used as variables to be solved for as part of the reconstruction. Explicit incorporation of these systematic deviations can make the reconstruction more robust to imperfections. With the parameters in hand, synthetic contrast-weighted images could be potentially generated by evaluating \autoref{eqn:mri:mxyt} with specific scan and sequence parameters, a technique known as Synthetic MR
\cite{bib:block2009,bib:tanenbaum2017kd}.

However, even though the individual terms are differentiable, the main downside is that the reconstruction problem is highly non-convex, thereby complicating the optimization.  The non-convexity results in increased computation time and dependence on initialization.
The problem requires a good initial guess of $\vec{\rho}$ and $\vec \theta$ in order to converge to a reasonable estimate. Another drawback is that mismatches between the model and the true acquisition can lead to error propagation in the estimated parameter maps. \revnum{2}{4}For example, radiologists have described the presence of flow artifacts, white-noise artifacts, and other artifacts in synthetic FLAIR contrast images \cite{bib:tanenbaum2017kd}, likely due to unmodeled effects such as flow. Partial voluming is also known to impair the estimation accuracy \cite{bib:block2009}. Errors in sampling trajectories due to eddy currents and gradient delays can also manifest as blurring and streaking artifacts.
These effects can potentially be reduced by expanding the signal model, e.g. to incorporate scanner non-idealities. The pulse sequence can also include navigator components to aid in the estimation.

\section{Relaxing the Multi-contrast Model}\label{sec:relax_mri}
Instead of fully incorporating MRI physics into the reconstruction, it is possible to use  relaxed constraints that are more amenable to optimization \cite{bib:doneva2010xi}. The natural dynamics of the MR signal are constrained by the Bloch equations, implying the existence of a low-dimensional manifold, which is non-linear in general.
As \autoref{fig:manifold_subspace_redo} illustrates, a linear subspace approximation to the manifold with higher dimensionality could provide a compromise between representation simplicity and size.
This may be attractive for a few reasons; namely, to maintain convexity and computational efficiency in the optimization, to decouple the reconstruction from the quantitative fitting, and to reduce propagation of model error. In addition, many applications do not require quantitative parameters, and instead rely on high-quality contrast-weighted images. This is the case for nearly all clinical diagnostic imaging. Unfortunately, contrast-weighted images derived from parameter maps are susceptible to error propagation due to unmodeled components, e.g. from partial voluming and flow effects, and have seen limited use clinically \cite{bib:tanenbaum2017kd}.
\myfig[manifold_subspace_redo]{The low-dimensional manifold representing the Bloch equations is captured by a linear subspace $\vec \Phi_K$ of larger dimension.}

\subsection{Subspace Constraint}
Continuing with the FSE sequence as a guiding example, \autoref{fig:fse_sim_white} shows the signal evolutions for an FSE simulation with a particular flip angle schedule. Despite differences in relaxation parameters, the signal evolutions for different tissues follow similar trends. This correlation implies low-dimensional structure; namely, the signal evolutions of different tissues form a low-dimensional subspace \cite{bib:liang2007ye,bib:huang2012oe, bib:zhao2014, bib:t2shmethod,bib:velikina2015jx}.
\myfigfull[fse_sim_white]{Forming a subspace based on signal dynamics, shown for FSE. (left) An ensemble of signal evolutions is drawn from a prior distribution. (middle) Due to correlation in signal dynamics, the data matrix is low-rank and (right) is well-approximated with PCA.}

Many approaches can be taken to design the subspace. Here we focus on simulating a set of training signals derived from EPG simulation. We assume a prior distribution $p(\vec \theta)$ is known and draw $L$ samples from the distribution to create the training signals. The prior distribution can be taken from known literature values, or from a conventional mapping procedure focused on a particular anatomy.

Consider a data matrix $\vec X \in\complex{T\times L}$ consisting of an ensemble of $L$ signal evolutions sampled at $T$ echo times. Each column in $\vec X$ represents
the signal evolution of a spin population with a particular $\vec \theta \in p(\bm \theta)$. Let $\vec \Phi\in \complex{T\times T}$ be an orthonormal
temporal basis, i.e.
\begin{align}
  \vec X = \vec\Phi\vec\Phi^H\vec X.
  \label{eqn:t2sh:span}
\end{align}
The goal is to design $\vec\Phi = \begin{bmatrix}\bm\varphi_1 & \cdots & \bm\varphi_T\end{bmatrix}$ and a $K-$dimensional subspace, \sloppy
$\vec\Phi_K = \spn\{\vec\varphi_1, \dots, \bm\varphi_K\}$, such that
\begin{align} \norm{\vec X - \vec \Phi_K\vec\Phi_K^H \vec X} < \epsilon,
  \label{eqn:t2sh:span2}
\end{align}
where $\epsilon$ is a modeling error tolerance.
The choice of norm in \autoref{eqn:t2sh:span2} will affect the chosen subspace and can be used to capture average, worst-case, and other error metrics.

When the Frobenius norm is used, the solution corresponds to the truncated singular value decomposition
\cite{bib:liang2007ye,bib:huang2012oe,bib:t2shmethod}, i.e. principal component analysis (PCA), where $\vec \Phi_K$ consists of the left singular vectors of $\vec X$ corresponding to the $K$ largest singular values. The principal component images, corresponding to the principal component vectors, are given by
\begin{align}
    \vec \alpha = \vec\Phi_K^H \vec x.
    \label{eqn:pca}
\end{align}
Based on \autoref{eqn:t2sh:span2}, we also have $\vec x \approx \vec \Phi_K \vec\alpha$.

\subsection{Reconstruction}
The subspace relationship can be incorporated into the reconstruction as additional prior knowledge in the form of regularization \cite{bib:velikina2015jx}:
\begin{align}
  \min_{\vec x} \frac{1}{2}\norm{\vec y - \vec E\vec x}_2^2 +
  \frac{\mu}{2}\norm{\vec x - \vec\Phi_K\vec\Phi_K^H \vec x}_2^2
  + \lambda \mathcal{R}\left(\vec x\right),
  \label{eqn:mocco}
\end{align}
where $\mu>0$ controls the degree of subspace modeling. Alternatively, the subspace can be used as an explicit constraint, i.e.,
\begin{nalign}
    & \underset{\vec x}{\text{minimize}}
& & \frac{1}{2} \norm{\vec y - \vec{E} \vec{x}}_2^2
+ \lambda \mathcal{R}\left(\vec x\right)\\
& \text{subject to}
& & \vec x = \vec \Phi_K \vec \alpha.
\label{eqn:subspace_constrained_recon}
\end{nalign}
Although the analysis form described by \autoref{eqn:mocco} is a more faithful representation, the reconstruction still requires solving for $TN$ parameters in $\vec x$, in addition to introducing a hyper-parameter $\mu$. In contrast, the synthesis form used in \autoref{eqn:subspace_constrained_recon} introduces explicit model error but is significantly more computationally efficient.

When using the hard constraint, we can now solve for the subspace coefficient images directly \cite{bib:t2shmethod},
\begin{align}
  \min_{\bm\alpha} \frac{1}{2}\norm{\vec y - \vec E\vec\Phi_K\bm\alpha}_2^2 + \lambda\mathcal{R}\left(\vec \alpha\right),
  \label{eqn:t2shrecon}
\end{align}
followed by back-projection: $\hat{\vec x} = \vec\Phi_K\vec\alpha$. This is a significant dimensionality reduction! Instead of solving for $TN$ variables, we only need to solve for $KN$ variables. In addition, when the normal equations are used in an iterative optimization (as is the case for many first-order iterative algorithms), we can take advantage of the commutativity of the subspace operator and the Fourier transform, as the former only operates on the parametric dimension and the latter only operates on the spatial coordinates \cite{bib:t2shmethod}:
\begin{align}
  \vec P \vec F \vec \Phi_K & = \vec P \vec \Phi_K \vec F\\
  \implies  \vec \Phi_K^H \vec F^H \vec P \vec F \vec \Phi_K &= \vec F^H\vec \Phi_K^H\vec P\vec \Phi_K\vec F \\
  &= \vec F^H \vec \Psi_K \vec F,
  \label{eqn:stkern}
\end{align}
where $\vec\Psi_K \in \complex{KM\times KM}$ is a block-wise diagonal operator with a $K\times K$ symmetric block  for each spatial frequency point. In other words, the optimization does not have to perform any computation in the ambient space, and the complexity grows with $K$, independent of $T$.

Compared to model-based quantitative mapping, the subspace-constrained forward model is convex and easier to solve. In addition, the subspace is less sensitive to model error, as it does not strictly impose a specific physical model. As an example, a voxel with partial voluming will contain a linear combination of signal evolutions, $\vec x(\vec r) = a_1\vec x^{(1)}(\vec r) + a_2\vec x^{(2)}(\vec r)$, comprising different tissue parameters $\vec \theta^{(1)}$ and $\vec \theta^{(2)}$, which is inconsistent with \autoref{eqn:mri:mxyt}. In contrast, if $\vec x^{(1)}$ and $\vec x^{(2)}$ are separately represented by the subspace, then so is their combination.

The main drawback to the subspace formulation is that the problem is not reduced to its intrinsic dimension governed by $\vec \rho$ and $\vec \theta$. However, in practice a subspace size of $K<5$ is practical for applications even when the ambient dimension is in the hundreds \cite{bib:t2shmethod,bib:zhao2017}.
A second drawback is that the subspace can represent points off the manifold, that are not physically meaningful. Thus, data inconsistencies can manifest as inaccurate images after back-projections.

\subsection{Quantitative Mapping}
A straightforward approach to quantitative mapping from reconstructed multi-contrast images is to perform a voxel-wise non-linear least squares fit based on the signal model. Given the reconstructed image $\hat{\vec x}$ at voxel $\vec r$, this amounts to solving
\begin{align}
  \underset{\rho(\vec r),\bm\theta(\vec r)}{\arg\min}\frac{1}{2}\norm{\hat{\vec x}(\vec r) - \rho(\vec r)\vec f(\bm \theta(\vec r), \vec u(\vec r))}_2^2.
  \label{eqn:nnls}
\end{align}

As with \autoref{sec:model_based_mri}, it is also possible to solve for RF field inhomogeneity by additionally solving for the non-negative scalar $\eta$ that multiplies the flip angles: $\vec u \rightarrow \eta\vec u$.
This formulation also covers dictionary-based methods which have grown in recent popularity with the advent of MR Fingerprinting \cite{bib:ma2013vg}. Rather than solving a continuous non-linear least squares problem, dictionary-based fitting utilizes a grid search across the parameter space and is equivalent to matched filtering. Optimization-based methods can be used as well.%

\revnum{3}{5}When the reconstruction strictly enforces a subspace, there is a known model error between the reconstruction $\hat{\vec x}$ and the signal evolution $\vec{f}$, and we do not expect them to match even in the absence of noise and under-sampling artifacts.
To improve parameter estimation when using a subspace-constrained reconstruction, we can solve for the parameters directly in the subspace:
\begin{align}
\underset{\rho,\vec \theta}{\arg\min}\frac{1}{2}\norm{\hat{\bm \alpha}(\vec r) - \vec \Phi_K^H \rho\vec f(\bm \theta, \vec u(\vec r))}_2^2.
  \label{eqn:nnls_subspace}
\end{align}

Since $\vec f$ is differentiable with respect to $\vec \theta$ and $\vec u$, first-order and second-order solvers can be used. \autoref{fig:t2map} shows an example of solving \autoref{eqn:nnls_subspace} following a $T_2$\xspace Shuffling reconstruction\footnote{Example code available on \url{https://eecs.berkeley.edu/~mlustig/Software.html}} \cite{bib:t2shmethod}, where $K=3$ and  $\bm \theta = (|\rho|, \angle \rho, T_2)$. The non-linear least squares was solved using the Trust Region Reflective algorithm included in the Python Scipy package, and the Jacobians were calculated using the adjoint states method \cite{bib:sbrizzi}.
\revnum{1}{7}
\myfig[t2map]{Example showing voxel-wise parameter fitting following a $T_2$\xspace Shuffling reconstruction using non-linear least squares directly in the subspace. The principal component images (magnitude and phase) were first reconstructed with the convex formulation \autoref{eqn:t2shrecon}, and the complex values were fit to the physical model using the EPG formalism.}

\revnum{2}{6}When solving under-determined inverse problems, it is important to recognize the bias introduced through modeling assumptions and regularization. In particular, the optimization problem \autoref{eqn:subspace_constrained_recon} introduces two forms of bias: model error due to the subspace constraint, and error due to the regularization. The subspace constraint leads to a straightforward tradeoff between model error and noise amplification: noise standard deviation increases with $\sqrt{K}$, where $K$ is the subspace size \cite{bib:t2shmethod}. A small subspace will also reduce sensitivity in the parameter space, and manifests as reduced contrast and parameter-dependent noise amplification, as all voxels are pushed to the same curve fit \cite{bib:huang2012oe}.
Bias due to regularization can also be an issue in reconstruction and parameter fitting, but its impact depends on the specific regularization used. With non-linear regularization such as $\ell_1$ commonly used in CS, the tradeoff can be harder to quantify. Blocking and smoothing artifacts are common when using wavelet regularization and total variation due to loss of high-frequency content \cite{bib:lustig2008pe}. Low-rank regularization and can lead to washed-out contrast due to loss of spectral information, and its multi-scale variants can lead to blocking artifacts due to their translation-variant structure \cite{bib:frankmlr}. In general, under-regularization can cause residual incoherent under-sampling artifacts and noise amplification, while over-regularization can lead to blurring \cite{bib:block2009,bib:doneva2010xi,bib:huang2012oe,bib:zhao2014,bib:zhao2014b,bib:velikina2015jx,bib:feng2016zk,bib:t2shmethod,bib:zhao2017,bib:jiang2018pe}.

\section{Choosing Scan Parameters}\label{sec:scan_params}
In addition to incorporating MRI physics in the reconstruction, using knowledge of the physics to optimize the MRI pulse sequence could potentially improve the quality of the acquired raw data. Designing an MRI pulse sequence consists of two main components: (\textit{i}) designing the pulse waveform and timing to guide the signal evolution, and (\textit{ii}) designing the spatial encoding to appropriately sample k-space. Initial work focused on designing sequences with simple goals in mind such as maximizing contrast and limiting blurring. A recent trend in computational MRI is to use optimal control to motivate scan parameter selection \cite{bib:sbrizzi}, and incoherent sampling properties to motivate k-space sampling \cite{bib:lustig2008pe}.
Optimal experiment design has also been applied to spatial encoding under a synthesis sparsity model \cite{bib:haldar2018nj}.
Importantly, both spatial encoding and sequence parameters should be designed with a specific reconstruction in mind. %

\subsection{Optimizing Sequence Parameters}
Many approaches to optimizing sequence parameter selection involve the use of a Cramer-Rao lower bound (CRLB), which imposes a lower bound on the variance of an unbiased estimator. To simplify the optimization, the spatial encoding is decoupled from the sequence evolution, and a 1D experiment is used for optimization:
\begin{align}
    \vec y = \vec f\left(\vec \theta, \vec u\right) + \vec w.
\end{align}
Then, the Fisher information matrix, $\vec I(\vec \theta;\bm u)$, which measures the sensitivity of $\vec \theta$ captured by $\vec y$ given sequence parameters $\vec u$,
can be used to lower bound the variance of an unbiased estimator.
The sequence parameters can then be optimized to minimize this lower bound on the variance of the biophysical tissue parameter estimates.
This approach has been used in \cite{bib:Nataraj2016} to develop a general framework to optimize sequence parameters from combinations of pulse sequences for precise estimation of $T_1$\xspace and $T_2$\xspace parameters jointly. It has also been used to determine sequence parameters for MR Fingerprinting to obtain maximal SNR while not violating physics-based MRI constraints \cite{zhao2016optimal}.

Another optimal control design based approach uses the EPG formalism to develop a model for the signal evolution in terms of various sequence and tissue parameters \cite{bib:sbrizzi}. An optimization problem is then formulated according to the goal, e.g.\ maximizing signal intensity or minimizing RF power, and solved using the adjoint states method. This approach can flexibly incorporate other objective functions, e.g. to minimize the CRLB, similar to methods discussed above.
\autoref{fig:t2shcrlb} shows an example of optimizing the RF flip angles in an FSE experiment given RF power constraints \cite{bib:keerthivasan2019}. The objective function maximized the component in the Fisher Information matrix corresponding to the variance with respect to $T_2$\xspace estimation. Constant flip angles that achieve the same RF power limit are shown in comparison.
\myfigfull[t2shcrlb]{Comparison of a constant flip angle schedule vs. optimized flip angles for an FSE experiment given a maximum RF power constraint (left).  The CRLB was optimized with respect to $T_1=1000$ ms and $T_2=100$ ms, and the received signal was higher for the optimized flip angles (middle). Even though the optimization targeted a single relaxation value, the bound was lower across a uniform range of $T_2$\xspace values (right).}

Data-driven methods for sequence parameter optimization are also emerging. With these methods, a certain set of sequence parameters are chosen, the resulting time-evolution of the signal for these parameters is either simulated or determined experimentally, and finally an image is reconstructed. The sequence parameters are then updated according to the reconstruction error. Compared to the optimal experiment based approaches, these methods have the advantage of incorporating spatial encoding and reconstruction into the evaluation.

\revnum{3}{6}When choosing sequence parameters, it is important to understand the impact of system imperfections on the optimized sequence. These imperfections can often be incorporated into and mitigated in the optimization framework, though they may greatly increase the computational complexity. For example, to reduce sensitivity to RF transmit field inhomogeneity, one may optimize the average (or min-max) CRLB across spins experiencing a distribution of RF transmit field inhomogeneities by simulating each spin independently. This quickly increases in complexity as additional components are considered.

\subsection{Optimizing Sampling Pattern}
Conventional methods to determine CS sampling trajectories are based on exploiting the incoherence conditions for CS by the use of a transform point spread function (TPSF), which determines the leakage of one transform coefficient to another transform coefficient due to subsampling \cite{bib:lustig2008pe}.
A sampling pattern can be determined by a Monte-Carlo design procedure where a variable density pattern is obtained by randomly drawing indices using a probability density function, the TPSF of the obtained trajectory is calculated, and the procedure is repeated, choosing the pattern with the lowest peak interference.
The reconstruction is performed according to \autoref{eqn:lasso} with a sparsity-promoting regularizer on the transform coefficients. Thus, incoherence in the transform domain is a good metric for sampling pattern quality.

Several data-driven approaches have been developed in order to improve on the current incoherence-based methods. These data-driven approaches use a learning-based framework to find the best sampling pattern for a set of training signals. In particular, methods that optimize a sampling pattern for a specific reconstruction rule and anatomy \cite{bib:haldar2018nj} as well as methods that jointly optimize a sampling pattern and a reconstruction strategy have been developed \cite{bahadir2019learning}. In this way, physics-based constraints in the reconstruction are implicitly propagated through to the sampling pattern optimization.\revnum{1}{8}

Although various control design based and data-driven strategies have been developed to optimize one of scan parameters, sampling trajectories and reconstructions, the ultimate goal of developing methods that optimize parameter, trajectory and reconstruction in conjunction remain an area for future work.

\section{Summary}
Since the advent of MRI, physics-based knowledge has been used to recover and understand the signal dynamics and tissue parameters. Many early works have leveraged physical constraints to mitigate image artifacts due to systematic errors and derive quantitative maps. However, long scan time and the lack of sophisticated reconstruction algorithms have prevented the clinical adoption of these techniques.

In the last decade, compressed sensing and other computational imaging approaches have transformed the landscape of what is possible with MRI. Scan time can be appreciably reduced by leveraging natural image statistics in the reconstruction. Many advanced reconstruction algorithms have been developed for this express purpose. Using these numerical tools and combining sparsity-based modeling, it is now feasible to run physics-constrained computational MRI methods in the clinic with reasonable scan and reconstruction times \cite{tamir_targeted_2019}.
By incorporating the physical dynamics due to both tissue-specific and scanner-specific parameters, acquisitions and reconstructions can be designed in tandem to work across a broad patient population in a robust manner.\revnum{2}{2}

The methods introduced in this article provide a framework for modeling the MR dynamics, modifying the acquisition to account for the signal evolutions, and incorporating them into the reconstruction. Although we limited our focus to only modeling dynamics due to relaxation effects and RF pulses, there are many other contrast mechanisms and scanner controls that can be accounted for. On the tissue characterization side, these include diffusion, water-fat imaging, susceptibility, spectroscopy, pharmacokinetics, magnetization transfer, chemical exchange, and blood flow. On the imaging system side, these include field inhomogeneity, eddy current effects, gradient delays, temperature, and more. Each of these components can be modeled jointly, but will greatly increase the dimensionality of the problem. Model-based quantitative imaging reduces the inverse problem to recovering the intrinsic parameters, but careful modeling and simulation must be used to avoid artifacts due to model error. In contrast, subspace constraints and other low-dimensional representations can be used to flexibly capture the dynamics without making strong assumptions.

Many exciting modeling techniques are emerging in the signal processing community, and are great candidates at improving physical modeling in MRI. On the other hand, purely data-driven approaches have grown in popularity, in the hopes of learning the signal characteristics from real data.
\revnum{2}{7}With the growing trend of applying deep learning to inverse problems, there is great promise to incorporating additional physics-based constraints directly in the learning in order to restrict the feasible solution space, reduce the dependence on large training data sets, and model effects not described by the simplified physics \cite{bib:modl}. Like compressed sensing, the rapid empirical progress in deep learning-based imaging should be suitably counter-balanced with theoretical guarantees to guide its use in clinical settings.
Combining physical modeling with data-driven learning is an active area of research in the signal processing community in general, and in the MRI field in particular.

\section*{Acknowledgements}
This work was supported by NIH grant R01EB009690, Sloan Research Fellowship, Bakar Fellowship, and research support from GE Healthcare. The authors thank Volkert Roeloffs, Zhiyong Zhang, and Karthik Gopalan for valuable comments.

\section{MRI Background}
\subsection{MRI Physics}
The MRI system can be approximated by a discrete linear dynamical system based on the Bloch equations, shown in \autoref{fig:mri_dynamical_system}. An aggregate of spins in each 3D spatial position (called a voxel) creates a net magnetization that is initially aligned with MRI's strong magnetic field, and evolves based on intrinsic biophysical tissue parameters and user control inputs consisting of radio-frequency (RF) and magnetic field pulses. RF pulses act to rotate or \textit{excite} the magnetization vector away from the main field direction, effectively mixing its longitudinal and transverse components.
The strength and duration of the RF pulse determines the degree of rotation, or \textit{flip angle}, experienced. The transverse component of the magnetization is sensed by nearby receive coils through Faraday's Law.
\myfig[mri_dynamical_system]{The MRI dynamical system. An aggregate of spins at each voxel position creates a net magnetization that evolves based on user pulse sequence control inputs and intrinsic tissue parameters.}

When no control inputs are played, the magnetization will return to the equilibrium state at a rate governed by the tissue's relaxation parameters. The longitudinal component will exponentially recover to its initial state with time constant $T_1$\xspace, and the transverse component will exponentially decay to zero with time constant $T_2$\xspace.
As \autoref{fig:Bloch_Toilet} shows, relaxation can be understood through analogy with a toilet, where the water in the tank represents longitudinal magnetization, water in the bowl represents transverse magnetization, and the toilet flush is an RF waveform generator. When the toilet is flushed (excitation), water transfers from the tank to the bowl, producing a detectable signal. As the toilet bowl drains ($T_2$\xspace relaxation), the tank refills with water ($T_1$\xspace recovery). Successive flushes will transfer new water from the tank to the bowl.
Neglecting other system effects, the magnetization is given by
\begin{align}
  M_z(t) &= M_0\left( 1 - e^{-\frac{t}{T_1}}\right), \label{eqn:t1} \\
  M_{xy}(t) &=  M_{0}e^{-\frac{t}{T_2}} \label{eqn:t2},
\end{align}
where $M_{z}(t)$ and $M_{xy}(t)$ are the longitudinal and transverse magnetization components at time $t$ after the excitation, respectively, and $M_{0}$ is the magnetization at thermal equilibrium.

\myfig[Bloch_Toilet]{Signal relaxation visualized through the toilet analogy. The longitudinal and transverse time constants, referred to as $T_1$\xspace and $T_2$\xspace, respectively, are intrinsic physical parameters and differ across tissue types.}

In general, the transverse magnetization can be described by
\begin{align}
  M_{xy}(\vec r, t) = \rho(\vec r)f_t(\vec{r}),
  \label{eqn:mri:mxyt}
\end{align}
where $\rho(\vec r)$ is the intrinsic amount of magnetization at position $\vec r$ (called proton-density, or PD), and $f_t(\vec r)$ is a spatio-temporal signal evolution that depends on the tissue and pulse sequence parameters.
By manipulating the sequence controls and sampling the magnetization at appropriate times, many different types of image contrasts can be captured. \autoref{fig:brain_multicontrast} shows four common image contrasts based on PD, $T_1$\xspace, and $T_2$\xspace.
\myfig[brain_multicontrast]{Image contrasts based on PD, $T_1$\xspace, and $T_2$\xspace produced by manipulating pulse sequence timing.}

\mysubsection{Frequency Domain Formulation}
The MRI signal samples spatial frequencies of the transverse magnetization, and the linear relationship is described by the integral
\begin{align}
  s(t) = \int_{\vec r} M_{xy}(\vec r, t) e^{-j 2 \pi \vec{k}(t)^\top\vec r} d\vec{r} + w(t),
  \label{eqn:mri:signal_equation}
\end{align}
where $s(t) \in \mathbb{C}$ is the acquired signal at time $t$ after the RF excitation, $M_{xy}(\vec{r}, t)$ is the transverse component of the magnetization that depends both on spatial position $\vec r$ and time $t$, $\vec{k}(t)$ is a trajectory of 3D frequency points and $w(t)$ is complex-valued white Gaussian noise. The symbol $(\cdot)^\top$ is used to denote the transpose operation.
Because the spatial frequency points are denoted as $\vec k$, MRI acquisition is often described as sampling in \textit{k-space}.
After change of variables and discretization, the MR signal equation is represented as the linear operator
\begin{align}
  \vec y = \vec{E x} + \vec w,
  \label{eqn:forward_model}
\end{align}
where $\vec x \in \complex{N}$ is the image, $\vec y \in \complex{M}$ are the acquired data in k-space, $\vec E \in \complex{M\times N}$ represents the spatial encoding, and $\vec w \in \complex {M}$ is noise.
If each measurement is acquired at time $t=\mathrm{TE}$ (called the \textit{echo time}) after RF excitation, temporal magnetization effects can be ignored and the forward model is
\begin{align}
  \vec E = \vec{PF},
  \label{eqn:forward_model_simple}
\end{align}
where $\vec F$ is a Fourier transform operator and $\vec P$ is the operator that selects the acquired k-space samples.
\autoref{fig:forward_model} illustrates this forward model relationship. In the case of 3D Cartesian imaging, the acquired samples are represented as points on a 2D grid, though more complicated, non-Cartesian sampling trajectories are also possible \cite{NONCART}.

\myfig[forward_model]{The discretized MRI forward model given by $\vec y = \vec{PFx}$. The MR image is Fourier transformed and measurements are acquired in k-space. In the case of 3D Cartesian imaging, the acquired samples are represented as points on a 2D grid.}

It is common to receive the MRI signal in parallel from multiple receive coils placed around the object. In particular, each coil element is sensitive to a local region of the underlying object, and can be represented as the underlying image multiplied with spatially varying coil sensitivity maps. Concretely, with multi-channel coil arrays, $\vec y \in \complex{CM}$ consists of multiple measurement vectors
\begin{align}
  \vec y_j = \vec{PFS}_j \vec{x} + \vec w_j, \quad j=1,\dots,C
  \label{eqn:forward_model_multicoil}
\end{align}
where $\vec S_j \in \complex{N\times N}$ is a diagonal matrix with the spatial sensitivity values on the diagonal, and $\vec w_j \in \complex{M}$ is the noise at the $\nth{j}$ coil, respectively, and $C$ is the number of coil elements in the array. Because the received signals are acquired in parallel, the sampling operator $\vec P$ is shared among all the channels. Parallel imaging uses the multi-channel measurements to under-sample k-space while maintaining a well-posed reconstruction.

\subsection{Signal Dynamics}
For many applications, it is not practical to wait for relaxation effects to ``reset'' the spin system before acquiring each measurement, as the effective wait time or \textit{repetition time} (TR) can be multiple seconds. This would translate to hundreds of hours of scan time. To overcome this limitation, several samples can be collected in a single TR to accelerate the acquisition. This is the basis of the RARE technique \cite{bib:hennig1986}, also known as fast spin-echo (FSE) or turbo spin-echo (TSE). However, the transverse magnetization evolves while sampling, breaking the simple Fourier transform relationship.
When left unmodeled, the signal dynamics modulate the acquired data and cause an undesirable filtering effect on the image. In the case of FSE, this leads to spatially varying image blurring due to $T_2$\xspace relaxation, shown in \autoref{fig:fse_blurring}.
\myfig[fse_blurring]{Effects of collecting samples at multiple TEs with FSE. (left column) A 90$^\circ$ excitation pulse (red rectangle) followed by 180$^\circ$ refocusing pulses (gray rectangles) will lead to exponential decay due to $T_2$\xspace relaxation. (middle column) When a single sample is acquired in each TR duration (red circles), a sharp image can be reconstructed. To accelerate the acquisition, measurements can be acquired between each refocusing pulse, shown as digital-to-analog (DAC) blocks. However, k-space will be modulated by $T_2$\xspace relaxation based on its acquisition time (multi-colored circles), leading to image blurring (right column). Since $T_2$\xspace depends on position, each voxel will exhibit different amounts of blurring.}

The linear forward model can be extended to represent the signal dynamics in the FSE experiment by including additional image states, $\vec x = \begin{bmatrix} \vec x_1 & \cdots & \vec x_T \end{bmatrix}$, where $\vec x_i \in \complex{N}$ is the image at the $\nth{i}$ imaging time and $T$ is number of time points. The forward model, illustrated in \autoref{fig:temporal_forward_model}, includes different encoding operators for time point $i$ and coil $j$:
\begin{align}
  \vec E_{ij} = \vec P_i \vec F \vec S_j.
  \label{eqn:temporal_forward_model}
\end{align}
The same concept can be used to include additional dimensions such as motion states due to cardiac or respiratory motion \cite{bib:lustig2007hl,bib:jiang2018pe}, spectral states for MR spectroscopy, dynamic contrast enhancement states \cite{bib:feng2014sm,bib:zhang2015iq}, and many others.
  \myfig[temporal_forward_model]{The forward model extended to relaxation dynamics for FSE. Each image represents an sampling time during the FSE acquisition, $T=4$ time points in this example. The forward model includes the sensitivity map operator ($C=2$ channels), the Fourier transform, and sampling. Each time point can be sampled differently, for instance using a center-out ordering of k-space.}

In the following sections, we see how modeling the signal dynamics allows us to reduce these artifacts as well as solve for the intrinsic tissue parameters that govern the dynamics.

\section{Pulse Sequence Design}\label{subsec:system_design}
Designing an MRI pulse sequence consists of two main components: (\textit{i}) designing the pulse waveform and timing to guide the signal evolution, and (\textit{ii}) designing the spatial encoding to appropriately sample k-space. Initial work focused on designing sequences with simple goals in mind such as maximizing contrast and limiting blurring \cite{bib:mugler2014gc}. A recent trend in computational MRI approaches is to use optimal experiment design to motivate scan parameter selection \cite{TODO}, and incoherent sampling properties to motivate k-space sampling \cite{bib:lustig2008pe}. Importantly, both should be designed with a specific reconstruction in mind. Although optimal experiment design in principle can be applied to spatial encoding as well, it leads to a prohibitively expensive optimization and remains an unexplored area of study. More recently, data-driven approaches that optimize (\textit{i}) and (\textit{ii}) jointly or separately have shown promise at achieving this goal.

\subsection{Sequence Parameter Selection}
The classical approach to creating signal contrast between different tissue types is to choose a specific imaging time that maximizes the difference in signal value between two target species. For example, focusing on $T_2$\xspace contrast given by \autoref{eqn:t2}, and defining the signal difference between two tissue species $A$ and $B$ as $M_{xy, A}(t) - M_{xy, B}(t)$, the imaging time that maximizes the difference is given by
\begin{align}
    t^\star = -\frac{\log\left(\frac{T_{2,A}}{T_{2,B}}\right)}{\frac{1}{T_{2,A}}-\frac{1}{T_{2,B}}}.
    \label{eqn:cnr}
\end{align}
The same concept is used to null tissue signal in inversion recovery sequences based on $T_1$\xspace (\autoref{fig:brain_multicontrast}).

Under the multiple measurement model \autoref{eqn:temporal_forward_model}, when multiple time points need to be chosen, an alternate formulation based on optimal experiment design is to maximize sensitivity to a particular parameter.

T1, T2 parameter maps are very useful biomarkers in monitoring the progression of various neurological and cardiac disorders.
Classical sequences such as inversion/saturation (IR/SR) or spin echo methods provide simple and reliable methods for estimating T1 and T2 maps, however, are limited by their long repetition time (TR).
On the other hand, steady-state sequences such as dual-echo steady-state (DESS) and Spoiled Gradient-Recalled Echo require a much shorter TR, greatly reducing scan time, and are sensitive to T1 and T2 variation. However, the T1, T2 cannot be directly extracted from the resulting signals requiring multiple acquisitions with different parameter sets to obtain T1, T2 maps.

Nataraj \textit{et al.} provide a general framework to optimize the sequence parameters (TR, TE and flip angle) of combinations of pulse sequences in order to provide accurate and precise parameter maps.
Inspired by the Cramer-Rao Bound (CRB), which bounds the minimum variance of an unbiased estimator, the set of scan parameters are optimized to minimize the worse-case maximum cost. The constructed optimization problem is solved with grid search.
The method is demonstrated on T1, T2 estimation from a 3 sets of different combinations of DESS and SPGR scans.
Numerical simulation results for optimized scan parameters of the 3 sets of combinations demonstrate that precise T1 and T2 parameters maps can be obtained.
Although the optimized scan profiles consisting of a combination of SPGR and DESS scans yield accurate results for the phantom dataset, these results do not hold up for the in vivo datasets, unless scan time is increased further.

Gradient-echo sequences can allow for faster acquisitions, however, after the preparation module signal decays gradually resulting in the initial contrast being reduced to a lower steady-state contrast level.
Two techniques to tackle this issue are segmented k-space acquisitions where a preparation module is repeated at the start of each segment, at the expense of longer scan times, and centric phase-encode ordering where low frequencies are acquired first after preparation, which results in reduced image resolution as higher spatial frequencies experience signal decay. This results in a tradeoff between image contrast and scan efficiency.

Cukur \textit{et al.} combine a post processing technique with centric phase-encode ordering and apply it to magnetization-prepared steady-state free precession (SSFP) angiography which improves contrast with no increase in scan time. The technique is based on the fact that in lower-extremity angiograms, the background mainly consists of skeletal muscle tissue with relatively lower spatial frequency content and the tissue of interest is thin blood vessels that have higher spatial-frequency content.
A compensation filter inversely proportional to the level of magnetization of the blood vessels is applied the acquired signal in the frequency domain, which compensates for the signal decay in blood vessels with no significant effect on the background muscle signal.
The compensation filter is combined with a sparsity-enforcing compressed sensing (CS) reconstruction, which can reduce the increased image noise due to the increased weighting of the higher spatial frequency components as a result of the compensation filter.
The highly compressible angiographic dataset is sampled with a variable-density random undersampling pattern in the k-space with decaying probability density towards higher spatial frequencies. The reconstruction is performed by solving a unconstrained nonlinear optimization problem that simultaneously enforces that the resulting reconstructed image is compressible and consistent with the acquired data.

Rapid Acquisition with Refocused Echoes (RARE) imaging uses an echo train consisting of echoes with changing phase encoding gradient allowing multiple lines in the k-space to be acquired in a single repetition time (TR). This sequence demonstrates that refocusing flip angles can be used to produce a higher signal intensity than the first echo as a result of constructive interference.
Busse \textit{et al.} present a technique to optimize the flip angles of the echo train to increase signal-to-noise ratio (SNR) of the resulting signal and produce a modulation transfer function (MTF) for different materials that improves resolution.
In addition they calculate a "contrast-equivalent" TE for each echo in the train, that matches the contrast that would have been achieved by a single echo, and show that the image contrast can be manipulated.
 The method proposed by the paper uses the inverse of the EPG algorithm ("Prospective EPG Algorithm") to find a sequence of flip angles that quickly establish psuedosteady-state, then keep the signal level at a constant rate for a certain number of echoes and finally start increasing the flip-angle linearly until the maximum flip angle is achieved.
In order to rapidly establish pseudosteady-state, a signal target for each echo in the echo train is determined that asymptotically approaches the signal target, $S_{target}$. This asymptotic technique can maintain a greater signal level and requires less RF power compared to alternatives including having a constant signal level starting from the first echo or using constant flip angle until maximum signal is achieved.
Other alternatives to increasing the flip angle linearly are considered including a Gaussian and half-exponential profile, however, the varying the angle linearly was able to best satisfy the resolution and SNR requirements. This technique shows that the flip angles of the echo train can be optimized to improve SNR and limit blurring artifacts due to relaxation of different materials.

Selective-inversion-recovery (SIR) is an imaging method that measures the quantitative magnetization transfer (qMT) or spin exchange between protons in free water and macromolecules by applying an RF pulse that selectively inverts the protons in free water.
The precision and accuracy of the qMT measurements from a SIR sequence depend on several parameters including the inversion recovery time ($t_i$) and the predelay time ($t_d$).
The signal measured by the SIR sequence can be expressed as a function of qMT parameters by adding the coupling between different pools to the Bloch equations.
Li \textit{et al.} introduce a CRLB based optimization technique that varies both $t_i$ and $t_d$ sampling points in order to maximize the qMT and T1 measurement precision efficiency separately for the worse case combination of parameters. Monte Carlo simulations are also performed by adding Gaussian noise to measure the uncertainties of the optimized parameters.
The results show that higher precision efficiencies can be obtained for both T1 and qMT measurements compared to previous methods and the uncertainties of the fitted parameters are reduced.

\subsection{Trajectory/Acquisition Strategies}

Various trajectory/acquisiton strategies have been explored to find the optimal k-space trajectory for a given task. Halder \textit{et al.} formulate an exponential decay based parametric contrast model of an image sequence consisting of $L$ different frames acquired with different image acquisition parameters. Each frame is thus a different contrast modulation of the same underlying image image. They examine the best trajectory in order to achieve high-resolution image frames as well as parameter maps. Both a CRB-based theoretical analysis as well as phantom experiments are performed to compare several 2D Cartesian sampling schemes and determine which scheme yields the most accurate reconstruction.

Lustig \textit{et al.} determine a more efficient k-space trajectory by employing the fact that the MR image has transform sparsity, a sparse representation in a fixed mathematical transform domain.
If an image exhibits transform sparsity, and if undersampling in the frequency space results in incoherent artifacts in the image space, then the image can be recovered from the undersampled frequency space data with a nonlinear recovery method.
The sparsity of MR images are displayed by performing a sparsifying transform on fully sampled images and reconstructing the images from various subsets of the largest transform coefficients.
Transform Point Spread Function (TPSF), which determines the contribution of one transform coefficient to another transform coefficient, was used to measure the incoherence of the sampling trajectory.
A sampling pattern can be determined by a Monte-Carlo design procedure where you determine a variable density pattern obtained by randomly drawing indices using a probability density function, calculate the TPSF of the obtained trajectory and repeat the procedure several times. Choose the pattern with the lowest peak interference.
The reconstruction is performed with a constrained optimization problem that involves maximizing the sparsity of the reconstructed image (using the l1 norm of the sparsifying transform of the image) while having a constraint that enforces data consistency. \textbf{add equations here} An additional Total Variation term can be added to the objective which enforces the sum of absolute variations in the image to be smaller. This allows the image to be sparse in the transform space and by finite differences simultaneously.

Building on (cite SPARSE MRI), OEDIPUS is a principled method for designing sampling patterns for sparsity-constrained MRI which accounts for coil geometry, anatomy, image contrast and can be reused for different subjects.
OEDIPUS uses a design based on the CRB bound on the covariance of an unbiased estimator. The CRB-based design for a unconstrained linear model, $d=Af+n$ where $f \in \mathbb{C}^N$, is modified by informing the model that the image vector has a certain level of sparsity, i.e. $\norm{\Psi f}_{0}\leq S$ where $\Psi \in \mathbb{C}^{Q \times N}$ is a known sparsifying transform and $S$ is the sparsity level. The linear model can then be written as $d=A \Psi^{\dagger} U\widetilde{c}+n$. Here, $c=\Psi f \in \mathbb{C}^Q$ and $c=U\widetilde{c}$ where $\widetilde{c}$ contains the S nonzero entries of $c$ and $U$ is a matrix formed by concatenating the $S$ columns of a $Q \times Q$ identity matrix.
The CRB bound for the covariance matrix of the unbiased matrix in the linear model is

\begin{equation*}
    cov(\hat{f}(d)) \succeq (A^HA)^{-1}
\end{equation*}

The model, informed of the sparsity constraints becomes

\begin{equation*}
    cov(\hat{f}(d)) \succeq \Psi^{\dagger}U(U^H(\Psi^{\dagger})^HA^HA\Psi^{\dagger}U)^{-1}U^H(\Psi^{\dagger})^H
\end{equation*}

The addition of the sparsity constraint also leads to large computational complexity advantages since the expressions being inverted are only size $S \times S$.
CRB bound considered in the paper is only valid for unbiased estimators, however, the problem at hand is not necessarily unbiased.

\section{Physical modeling}\label{sec:physical_modeling}
\comment{Ke \\
This section will discuss low-dimensional modeling based on the physical properties of the signal. The section will focus on linear subspace design informed by the expected behavior of the signal given the system design in Section~\ref{subsec:system_design} \cite{bib:zhao2012,bib:huang2012,bib:huang2013,bib:velikina2013}.
}
As the MRI signal is generated by protons in the body under a strong static magnetic field $B_0$, the signal itself is constrained by nuclear magnetic resonance physics like radio frequency (RF) excitation and relaxation. In this regard, the physics properties of the signal can be leveraged for low-rank modeling and image reconstruction. Conventional imaging methods assume that the samples of the frequency domain $s(\textbf{k},t)$ has to satisfy Nyquist rate. However, it has been demonstrated that high-quality images can be reconstructed from undersampled k-space data with appropriate physics constraints. These constraints include partial separability and sparsity for dynamic imaging reconstruction\cite{bib:zhao2012}, low-rank model for relaxation estimation\cite{bib:huang2012,bib:t2shmethod}, the smoothness in the parametric dimension and the relevant reconstruction methods\cite{bib:velikina2013}.

On the other hand, by utilizing the physics model behind MRI using the Bloch Equation\cite{bib:bloch1946nuclear}, we can design powerful simulator tools for MR signal simulation, which can be used for optimizing scan parameters and sequence programming. Some useful simulator will be briefly discussed in this section. Additional modeling approaches, including union of subspaces \cite{bib:velikina2015reconstruction,Mandava:MagnResonMed:2018}, low-rank tensors \cite{bib:christodoulou2018} and manifold learning \cite{bib:nakarmi2017}, will be briefly introduced.
\subsection{Physics Model for Image Reconstruction}
For traditional imaging methods, Nyquist sampling rate in the frequency domain is required for image reconstruction. However, the signal evolution function of MRI follows physics model(bloch equation\cite{bib:bloch1946nuclear}) which can be used for a low-dimensional representation of the signal itself. With some physics based assumptions or
constrains, image reconstruction for dynamic imaging, quantitative mapping and multi-contrast imaging can achieve a higher undersamping rate. Among them, Zhao \textit{et al.} proposed a dynamic image reconstruction method from highly undersampled (\textbf{k},\textit{t})-space data to use partial separability (PS) and sparsity constraints jointly for real-time cardiac MRI\cite{bib:zhao2012}. The proposed model incorporates the complementary advantages of both PS and sparsity constrains into a single formulation by using sparsity constraint on PS regularization, which is shown as follow:
\begin{equation}
\hat{\vec{\rho}} = \arg \min_{\bm{\rho}\in \mathbb{C}^{NM \times 1}}\| \bm{d}-\bm{\Phi \rho}\|_2^2+R_r(\bm{\rho})+R_s(\bm{\rho}),
\end{equation}
where $\vec{d}\in\mathbb{C}^{D\times 1}$ contains the measured data, $\bm{\Phi}\in\mathbb{C}^{D\times NM}$ is the measurement matrix, consists of the spatial Fourier Transform and sparse sampling operator. $\bm{\rho}=\text{vec}(\bm{C})\in \mathbb{C}^{NM\times1}$ is the matrix and vector representation of the dynamic image $\rho(x,t)$. The penalty function $R_r(\cdot)$ and $R_s(\cdot)$ are used to incorporate the PS and sparsity constrains respectively.
Also, a globally convergent computational algorithm is described to solve the underlying optimization problem. By using these properties of the signal, the in-vivo experiments demonstrated the better performance than using partial separability and  sparsity constraints individually, and can be applied to more dynamic imaging applications.

Besides the dynamic imaging, physics model based reconstruction algorithm can also be incorporated to quantitative parameter mappings (T1, T2 mappings).Huang \textit{et al.} presented a model-based algorithm to reconstruct T2 mapping from highly undersampled data, by combining principal component analysis. The nonlinearity related to T2 estimation is remove by introducing PCA where the exponential decay can be accurately approximated by linear combination of a few pre-calculated principal components (PCs). To be more concrete, let $\overrightarrow{v}$ be a real vector representing a noiseless exponential $T_2$ decay curve for real $I_0$ values, the exponential decay model can be wriiten as:
\begin{equation}
    v_j=I_0\cdot \text{e}^{-TE_j/T_2}.
\end{equation}
Let $L$ be the number of principle components, the $T_2$ decay curve can be approximated by the linear combination:
\begin{equation}
    \overrightarrow{v} = \sum_{i=1}^{L}m_i\overrightarrow{b_i}+\overrightarrow{\epsilon},
\end{equation}
where $\overrightarrow{b_i}$ is the $i$th PC vector and $m_i$ is the corresponding weights.
Following this work, Huang \textit{et al.} further developed an algorithm for fast and accurate T2 estimation from highly undersampled multi-echo spin-echo data, leveraging a model-based reconstruction method based on slice-resolved extended phase graph (SEPG) model. The principal components (PCs) are calculated for given TE points using a training set of decay curves given by SEPG model. Also,Velikina \textit{et al.} proposed an accelerating MR parameter mapping using sparsity-promoting regularization in parametric dimension. The signal in the parametric dimension is described by a smooth function, which is viewed as a regularization term in a CS-type minimization problem. The images series obtained as the solution of this minimization problem is then fit to a parametric model to yield parameter maps.

\comment{Approaches to simulating the signal will be discussed, including Bloch Simulation \cite{bib:eliezer2014} and the Extended Phase Graph Algorithm \cite{bib:weigel2015}.}

\subsection{MR Signal Simulation tools}
In the past decades, approaches to simulate the MR signal has been well developed based on the physics models, including simulators based on Bloch equations\cite{bib:bloch1946nuclear} and the Extended Phase Graph\cite{bib:weigel2015}. As the fundamental physics model of NMR and MRI, the Bloch equations are a set of macroscopic equations used for calculating the nuclear magnetization as a function of time when relaxation times are present. By leveraging Bloch simulation, Ben-Eliezer \textit{et al.} presented a postprocessing approach to overcome the common penalties associated with multiecho protocols. Bloch simulations are repeated for a range of $\text{T}_2$ values and transmit field scales, yielding a dictionary of simulated echo-modulation curve (EMC), which is then used to determine the $\text{T}_2$ value by finding the most closely matches curve in the dictionary.

The extended phase graph is a powerful tool for simulating and understanding the magnetization response of MR sequences. The effect of gradients, radiofrequency (RF) pulses, relaxation, and motion phenomena during the MR sequence is characterized as the action of a few matrix operations on different magnetization states, which can be viewed as a useful tool for MRI sequences design.

\textbf{Figure 1} Show an example of how EPG works. Using the Fig 3 in the paper may be an option.

\subsection{Additional  Modeling  Approaches}
Besides the methods introduced in the previous section, additional modeling approaches, including union of subspaces \cite{bib:velikina2015reconstruction,Mandava:MagnResonMed:2018}, low-rank tensors \cite{bib:christodoulou2018} and manifold learning \cite{bib:nakarmi2017}, will be briefly introduced in this section. By leveraging the consistency in the temporal domain, Velikina \textit{et al.} introduced the model consistency condition (MOCCO) technique, utilizing temporal models to regularize reconstruction without constraining the solution to be low-rank\cite{bib:velikina2015reconstruction}.

\textbf{Figure 2} Show how mocco works, cite maybe fig 7 in the paper.

This method allows enforcing consistency with the data-driven pre-estimated model, which can improve the reconstruction of dynamic imaging. Following this method, Mandava \textit{et al.} presented an extension of the MOCCO framework to work with a union of local subspaces constraint for quantitative mapping. Magnetic resonance multitasking conceptualized different sources of image dynamics as different time dimensions, using a low-rank tensor (LRT) imaging method for cardiovascular imaging. Christodoulou \textit{et al.} developed a new LRT imaging strategy, where a factored tensor reconstruction approach that enforces the LRT in two efficient steps. Nakarmi \textit{et al.} proposed a novel kernel-based framework to learn and incorporate a nonlinear manifold into constrained reconstruction from reduced acquisition by using kernel principal component analysis (KPCA).

\section{Physics Constrained Image Reconstruction}
\comment{Frank
This section will describe a generalized inverse problem that captures both physics-based structure and image-based prior knowledge \cite{bib:csmri,bib:doneva2010,bib:zhao2014}. Different formulations will be described along with their tradeoffs, e.g.\ explicit constraints vs regularization, and their effects on physical parameter estimation.}

Given the acquired k-space measurements $\{\vec{y}_{ij}\}_{ij}$, there are in general two goals: one is to reconstruct the underlying image sequence $\{\vec{x}_j\}_j$, and the other is to recover the quantitative parameter maps $\vec{T}_1$ and $\vec{T}_2$. \TODO{rephrase}Certain applications, such as in musuclar skeletal imaging, desire the image sequence of different contrast for better visualization, while other applications, such as in cardiac imaging, wants the quantitative maps.  In this section, we will focus on image reconstruction.

Because of long acquisition time, it is desirable to shorten the number of samples. Such inverse problems are in general ill-posed and physical-based structures and image-based prior knowledge are used to help regularize the recovery problem.

Various works have been proposed and in general consider objective functions of the following form:
\begin{equation}
\begin{aligned}
  &\underset{\vec{x}}{\text{minimize}}~\frac{1}{2} \sum_{i = 1}^T \sum_{j} \| \vec{E}_{ij} \vec{x}_j - \vec{y}_{ij} \|_2^2 + \lambda \mathcal{R}(\vec{x})\\
  &\text{subject to}~\vec{x} \in \mathcal{C}
\end{aligned}
\end{equation}
where $\mathcal{C}$ is some physics-based constrained set, and $\mathcal{R}$ is some regularization function.

\myfig[physical_constraint]{Physics-based constraints}

One way of utilizing the physical forward model is to enforce the underlying image to lie on the exact physical constraint:
\begin{equation}
\begin{aligned}
  \mathcal{C} = \{[\vec{x}_1(T_1, T_2), \hdots, \vec{x}_T(T_1, T_2)] : (T_1, T_2) \in \mathbb{R}_+^2\}
\end{aligned}
\end{equation}

But there are a number of difficulties. First the set is non-linear and non-convex. Second, any systematic error in the forward model would cause huge errors in the result. Highly coherent. Hence, in practice, the . However for simple systems, there are some works. \TODO{Related works}.

Instead of computing the forward model explicitly, low-rank approaches can be used to reduce computation. Low-rank approaches rely on a partial separability model \cite{bib:zhao2012}: there is strong spatiotemporal correlation in the data between, for example, frames of a dynamic image. The model of $L^{th}$ order partial separability for dynamic imaging is:
\begin{equation}
s(k_x, k_y, k_z, t) = \sum_{l=1}^{L}{g_l(k_x,k_y,k_z)h_l(t)}
\end{equation}
where $h_l(t)$ are the temporal bases and $g_l$ are the spatial bases.
The acquired data can be stacked into a Casorati matrix that indexes over k-space and time. In undersampled acquisitions, this matrix will have missing entries. The goal is to fill in the missing entries of the Casorati matrix $C$ subject to $rank (C) \leq L$.
\textbf{Figure 2}: show the Casorati matrix (see googledoc).

This model becomes very useful for spatially and temporally incoherent acquisitions, such as Magnetic Resonance Fingerprinting.

Single evolutions can be simulated using the Bloch equation or EPG as described in Section~\TODO{section}. In particular, given an ensemble of spins sampled at $T$ echo times, $\vec{X} \in \mathbb{C}^{T \times N}$.

\begin{equation}
    [\vec{x}_1(T_1, T_2), \hdots, \vec{x}_T(T_1, T_2)]
\end{equation}

\begin{equation}
\begin{aligned}
  \mathcal{C} = \{\vec{x} = \vec{\Phi}_K \bm{\alpha}\}
\end{aligned}
\end{equation}

Some image based priors to incorporate:

Wavelet/gradient sparsity prior
\begin{equation}
\begin{aligned}
  \mathcal{R}(\vec{x}) = \| \vec{\Psi} \vec{x} \|_1
\end{aligned}
\end{equation}

Low rank prior
\begin{equation}
\begin{aligned}
  \mathcal{R}(\vec{x}) = \| \vec{x} \|_*
\end{aligned}
\end{equation}

\section{Quantitative Mapping}
\comment{Suma \\
This section will show how the reconstructed data can be used to fit quantitative parameters to the physical model, e.g.\ through direct estimation \cite{bib:block2009}, non-linear least-squares, dictionary-based matched filtering \cite{bib:zhao2017,bib:asslander2017}, and novel approaches that operate directly in the low-dimensional subspace defined by the physical model \cite{McGivney:IeeeTransMedImaging:2014}.
}

\subsection{Nonlinear least squares}
To estimate $T_1$ and $T_2$ maps, many approaches consider the reconstructed images and then fit quantitative parameters. Instead, it is possible to estimate $T_1$ and $T_2$ directly from k-space, without the intermediate step of reconstructing an image \cite{bib:block2009}. This is done via modeling the reconstruction, and solving iteratively for the parameters directly.

As an example, for reconstructing $T_2$ and $\rho$ (proton density) maps from an undersampled fast spin echo acquisition \cite{bib:block2009}, we would like to recover the pair $(\vec{\rho}, \vec{r})$ from the raw data $\vec{y}$, where $r = 1/T_2$ for each element in $\vec{r}$. We can model $(\vec{\rho}, \vec{r})$ as a sum of scaled delta functions, where each voxel in the image is scaled by $\rho_{i}$ and $r_{i}$. Let the functions that evaluate $\rho_{i}$ and $r_{i}$ at each point in space by denoted by $\varrho$ and $R$, i.e.:

\begin{equation}
	\varrho(\vec{x}) = \sum_i{\rho_i \cdot \delta(\vec{x} - \vec{x_i})}
\end{equation}

and
\begin{equation}
	R(\vec{x}) = \sum_i{r_i \cdot \delta(\vec{x} - \vec{x_i})}
\end{equation}

We can construct the following forward model of the $j^{th}$ acquired sample as:
\begin{equation}
    F_j(\vec{\rho}, \vec{r}, t,c) = \sum_{\vec{x} \in FOV}{ \varrho(\vec{x}) \cdot e^{-R(\vec{x}) \cdot t} \cdot C_c(\vec{x}) \cdot e^{-j\vec{x} \cdot \vec{k_j}}}
\end{equation}

where $C_c(\vec{x})$ is the coil sensitivity profile of the $c^{th}$ coil at position $x$, and $t$ indexes over echo times.

To solve for $(\vec{\rho},\vec{r})$, we can solve a nonlinear least squares objective. Generally, this is of the form:

\begin{equation}
    \theta \in \text{argmin}~\| \vec{E} \vec{x}_t (\theta) - \vec{y}_t \|_2^2
\end{equation}

In this case, the objective is:
\begin{equation}
	\Phi(\vec{\rho}, \vec{r}) = \frac{1}{2}\sum_t{ \sum_c{ \|F(\vec{\rho}, \vec{r}, t,c) - \vec{y_{t,c}} \|_2^2}}
\end{equation}

The objective can be minimized with the conjugate gradient algorithm. It is differentiable; however, it is highly nonconvex. Additionally, the solution may result in unrealistic map estimates because an area of k-space may be ill-conditioned. Therefore, regularization term $P(\vec{\rho}, \vec{r})$ can be added to ensure smoothness:
\begin{equation}
	\Phi(\vec{\rho}, \vec{r}) = \frac{1}{2}\sum_t{ \sum_c{ \|F(\vec{\rho}, \vec{r}, t,c) - \vec{y_{t,c}} \|_2^2}} + \lambda \cdot P(\vec{\rho}, \vec{r})
\end{equation}

\subsection{Dictionary matching}
Another approach is to generate a dictionary of possible signal evolutions and pattern-match the reconstructed image to it, voxel by voxel. This technique is commonly used in Magnetic Resonance Fingerprinting \cite{bib:ma2013}. In MRF, data are acquired with pseudo-random parameters such as flip angles and TRs. The trajectory is an undersampled variable density spiral so that there is both spatial and temporal incoherence, which leads to incoherent aliasing. The reconstruction estimates parameters by pattern-matching to a pre-computed dictionary, voxel by voxel. Since the signal evolutions are very similar to one another, it is possible to compress the dictionary via SVD \cite{McGivney:IeeeTransMedImaging:2014} or to solve for the spatial and temporal components directly via low-rank matrix factorization \cite{bib:zhao2017}. These approaches are not discussed in detail here, but are explained in a different article in this issue.

\subsection{Nonlinear least squares from reconstructed image}
\begin{equation}
    \theta \in \text{argmin}~\| \vec{x}(T_1, T_2) - \vec{x} \|_2^2
\end{equation}

\subsection{Dictionary matching}
\begin{equation}
    \theta \in \text{argmax}_\theta~\left< \vec{E} \vec{x}_t (\theta), \vec{y}_t \right>
\end{equation}

\subsection{Nonlinear least squares from k-space}
\begin{equation}
    \theta \in \text{argmin}~\| \vec{E} \vec{x}_t (\theta) - \vec{y}_t \|_2^2
\end{equation}

redpajama_set_name: RedPajamaArXiv


text:
Your Expat Community in San Francisco
Connect with fellow expats in San Francisco
Get information in our San Francisco guides
Exchange tips about expat life in San Francisco
Raul Gonzales
Thanks to InterNations, I quickly found useful tips of where to search for a flat. And I met great global minds at the same time.
Welcome to Our Expatriate Community of San Francisco
Hello and welcome to the InterNations Community of expats who left their hearts in San Francisco! Our community provides an invaluable library of comprehensive guides, advice, and insider tips especially created to support expats; you can find the answers to questions such as, "how can I access healthcare?", "when should I apply for a visa?", or "what is the climate like?". Day-to-day dilemmas are solved in our forum by experienced expats who have already negotiated the uncertainties of relocating to San Francisco. Our members are always happy to share their knowledge, from the best family-friendly festivals to finding a romantic restaurant in the Bay. InterNations opens the doors to America's gold rush boom town as you cross the Golden Gate Bridge to your new life.
With InterNations, you never stand alone. Whether you are moving to another country or going on a business trip somewhere — our communities around the world are there to make you feel at home wherever you are. Feel free to also check out our other vibrant InterNations Communities like the ones in Khobar (Eastern Prov.) or Dubai.
Events for Expats in San Francisco
San Francisco Expat Forums
More information on Expat Life in San Francisco Expat Forums
Connect with Like-Minded Expatriates in San Francisco
From Argentina, living in San Francisco
From Australia, living in San Francisco
Meet other expats like you in San Francisco
Join Our Exciting Events in San Francisco
Attend our monthly events and activities for San Francisco expatriates to get to know like-minded expatriates in real life.
InterNations SF Holiday Happy Hour @Minna Gallery
InterNations SF First Newcomers' Night of the Year!
Star Wars Rise of Skywalker - 155 min
An AI Paper and A Drink - Meta Learning-Holiday Edition!
Want to see all the upcoming events in San Francisco?
Exchange Tips about living in San Francisco
Get trustworthy advice and local insights from fellow members in our San Francisco expat forums.
SUR an Argentine affair
Jun 6, 06:52
Hola Diana, te puedo averiguar.Pero, si estas en Palo Alto te recomendaria Foothill College (Los Altos Hills) Cuando yo estuve ahi tenian buen equipo e instructores. Cualquier cosita escribime a …
Room in SF
Aug 19, 00:06
Hi Sonia, Great that you could manage, I booked myself in LQ in Berkeley. I am looking forward to a great week of work and fun. You have great time too. Cheers Rajesh
Needing Advice!
Hello Ashley, While I've never been an au pair, I think you have a pretty decent offer. Since your room/board, meals are free the $450 is yours to spend as you see fit. Here's the catch, San …
Got questions about living in San Francisco?
Salsa-Kizomba Group
Bay Explorers Group
Movie Night Out Group
Share your hobbies and interests with expats in San Francisco
Embark on a Busy New Life in San Francisco
Today, San Francisco, with a population of approximately 850,000 citizens, is among the most colorful American cities and expats can feel it at the moment they start walking around the streets. In keeping with its gold rush past, San Francisco is a center of finance with national and multinational banking companies stationed here. The city also houses an enormous technology industry covering research into biotechnology, aerospace, and internet development from companies such as Apple, Google, and Pinterest. Maybe because of its cosmopolitan spirit and tremendous opportunities for personal development, San Francisco is one of the most attractive towns for young expats. Possibly the greatest industry of San Francisco is tourism, but in a city of this size expats with the right experience will find a variety of paths to pursue. When looking for accommodation, you can choose between the colorful houses of the suburbs such as Walnut Creek, which is excellent for families, or choose to stay closer to the city in trendy loft apartments. For those looking for luxury, the immaculate Sea Cliff neighborhood provides sea views, access to beaches, the Golden Gate National Park and private swimming pools. And if you are worrying about your neighbors, stop now! People in San Francisco possess a very positive way of thinking and at the same time they smile and laugh a lot. So, you couldn't wish for better next-door neighbors. Speaking of traveling and public transportation, the city will not disappoint you, it has lots of steep hills, twisted streets with sharp turns and cable cars. In addition to vibrant nightlife and huge shopping malls, the city also hosts year-round colorful festivals such as jazz, Shakespeare and chocolate festivals, Oktoberfest by the Bay and the Chinatown Night Market Fair. Culture lovers can find Broadway musicals on tour in one of the many glamorous theaters, and sports fans can choose between basketball, baseball, football, hockey, and soccer. Let's not forget about all of the places and landmarks which this city offers to its citizens. If expats are looking for outdoor activities, then Ocean Beach, Baker Beach and Golden Gate Park are the perfect areas to spend your free time. Moreover, Coit Tower, which offers a 360-degree view of the city, is also an alternative for expats to get more information about San Francisco. Numerous museums, zoos, and theme parks mean that there's something for everybody, and don't forget to find an InterNations Group to share your new interests with. All in all, as one of the greatest journalists, Herb Caen, born and raised in San Francisco once said, "I hope I go to Heaven, and when I do, I am going to do what every San Franciscan does when he gets there. He looks around and says, 'It ain't bad, but it ain't San Francisco'."
InterNations is the largest global network for expatriates, with communities in 420 cities around the world. Meet fellow global minds in your city of residence or get to know local expatriates before you're moving or traveling abroad and never feel like a stranger. InterNations makes it easy for you to stay in touch with your expat friends, both online and offline, so you can grow your network of friends and contacts from all over the world. Whether they live in Jeddah or Shanghai, InterNations allows you to be close to them even if they are miles away.
Socializing at InterNations Events and Groups
At InterNations, we're committed to providing expats with valuable information. For us, the safety and security of every community is of great importance, so we value the key principles of honesty, ethical behavior and mutual respect. At our monthly social events, you can meet global minds from all over the world; take the opportunity to network or just relax and make new friends. Share your interests with others or discover something new with our groups dedicated to concerts, books, culture and softball, among others. With advice, support, and real life networking opportunities, it's reassuring to know that you're part of our international InterNations community.
Li Wang
As I already connected with several Chinese women before moving to San Francisco, I quickly felt at home here.
Other Communities in United States of America
Expats in Washington, D.C.
Expats in New York
Expats in Los Angeles
Expats in Chicago
Expats in Houston
Expats in Philadelphia
Expats in Boston
Expats in Miami
Expats in Detroit
Expats in Seattle
Expats in Atlanta
Expats in Dallas
Expats in Denver
Expats in Jersey City
Expats in Palm Beach
Expats in San Diego
Expats in Minneapolis
Expats in Las Vegas
Expats in Austin
Expats in Silicon Valley
Expats in Phoenix
Expats in Charlotte
Expats in Raleigh
Expats in Orlando
Expats in Portland
Expats in Pittsburgh
Expats in Omaha
Expats in Cincinnati
Expats in Columbus
Expats in Tampa Bay
Expats in New Orleans
Expats in Honolulu
Expats in Sacramento
Expats in Milwaukee
Expats in Salt Lake City
Expats in Norfolk
Expats in Orange County
Expats in Chattanooga
Expats in Spartanburg
Expats in Nashville
Like-Minded Expatriates in San Francisco
French in San Francisco
Germans in San Francisco
Indians in San Francisco
Brits in San Francisco
Brazilians in San Francisco
Italians in San Francisco
Australians in San Francisco
Russians in San Francisco
Canadians in San Francisco
Chinese in San Francisco
United States Guide Topics
Moving to San Francisco
Living in San Francisco
Working in San Francisco
Moving to the US
Relocating to the US
Visas & Work Permits in the US
Housing in the US
Healthcare in the US
Banks & Taxes in the US
Working in the US
Country facts in the US
Best Places to Live in the US
Moving from Canada to the US
Moving from Germany to the US
Forum Categories in San Francisco

redpajama_set_name: RedPajamaCommonCrawl


text:
As temperatures plummet again, Lansing could rewrite cold records
Region nears all-time mark of consecutive days below 20 degrees.
As temperatures plummet again, Lansing could rewrite cold records Region nears all-time mark of consecutive days below 20 degrees. Check out this story on lansingstatejournal.com: http://on.lsj.com/2lPMobL
Eric Lacy and Ken Palmer, Lansing State Journal Published 5:51 p.m. ET Jan. 3, 2018 | Updated 6:29 p.m. ET Jan. 3, 2018
Simply Money Advisors discuss how to save money on heating during the cold season.
Tony Triplett of Lansing and J.J., his 5-year-old nephew, took shortcuts Wednesday across an empty lot near Kalamazoo Street. Wednesday marked the ninth straight day in the region below 20 degrees.(Photo: Matthew Dae Smith/Lansing State Journal)Buy Photo
LANSING -- Think it's cold?
Prepare for the weather to get worse, and maybe make history.
As of Wednesday night, the Lansing area had experienced nine straight days of high temperatures below 20 degrees.
The area probably won't get above 7 degrees on Thursday or Friday, according to the National Weather Service's forecast.
Three tips to stay safe and warm on CATA
8 things you shouldn't leave in your car when it's cold
Dangerous wind chills expected in region
By Saturday, an anticipated mark of 12 straight days below 20 degrees would mark the second most in the region's history, said meteorologist Jared Maples, of the weather service's Grand Rapids office.
The region's single greatest streak, Maples said, is 14 days below 20 degrees set Jan. 16, 1912. The modern mark is 10, set on Dec. 24, 1989.
The week of Dec. 25 to Dec. 31 was the coldest on record in Lansing, Grand Rapids and Muskegon, the weather service said. During that seven-day stretch, temperatures averaged about 6 degrees in Lansing and Grand Rapids and 10 degrees in Muskegon.
Despite the bitter cold, officials haven't called for warming centers to open, as they did during the December 2013 ice storm, when tens of thousands of people were without power, some for as many as 11 days. That's because no one has seen a need to take that step, officials said.
"(If we have power), everyone kind of has heat," said Rob Dale, emergency planner for Ingham County. "As long as we've had power, we've never had any dramatic demands."
Dale said he, other area emergency management officials and representatives from various agencies — including homeless shelters, the Tri-County Office on Aging and CATA — periodically compare notes. So far, no one has asked for extra services, he said.
Joan Jackson Johnson, director of the Lansing's Human Relations and Community Services Department, said officials have not ruled out declaring a state of emergency if the need arises. That would prompt the opening of warming centers and other services, she said.
In Delta Township, the Lansing Mall and Mount Hope Church allow anyone to come inside to keep warm during normal business hours, said John Clark, the township's fire chief. The church hasn't reported anybody coming in for shelter from the cold, he said.
"Typically, you don't open warming centers on a large scale unless you have something else going on, like a power outage," he said.
Dale said people should check on their neighbors, particularly the elderly and folks who live alone, to make sure they are okay.
2017 gave us the coldest December 25-31 on record. Considering both the high and low temperatures, that 7-day stretch averaged about 6 degrees at Lansing and Grand Rapids, and 10 degrees at Muskegon. #wmiwx#miwx
— NWS Grand Rapids (@NWSGrandRapids) January 2, 2018
The day center at the Volunteers of America complex on Larch Street in Lansing has been "extremely busy," said Darin Estep, community engagement director at the VOA. And the overnight shelter has been pushed into overflow mode repeatedly, he said.
As of Tuesday, 21 nights had been declared cold-weather emergencies since Nov. 9, prompting the agency to lift normal limits and put out extra beds, Estep said.
"It's unusual for us to see this many overflow nights this early in the season," he said. "Here we are, just getting into January, and who knows what will happen from here."
Eric Lacy is a reporter for the Lansing State Journal. Contact him at 517-377-1206 or elacy@lsj.com. Follow him on Twitter @EricLacy. Contact Ken Palmer at (517) 377-1032 or kpalmer@lsj.com. Follow him on Twitter @KBPalm_lsj.
Record cold
Streaks of days below 20 degrees in Lansing area:
14 - set Jan. 16, 1912
10 - set Dec. 24, 1989
9 - set Jan. 3, 2018
9 - set Jan. 1, 15, 1999
9 - set Jan. 29, 1961
9 - set Feb. 13, 1899
Read or Share this story: http://on.lsj.com/2lPMobL
Police investigating woman's murder in south Lansing
Clifford, Michigan horse known for his art, has died
Scientists: Deadly heat to increase in Michigan
LGBTQ student group joins lawsuit over Williamston policies
State denies disaster aid for Portland flooding


Double decker bike rack is ideal for places like supermarkets, apartments, schools etc for it saves spaces and securely protect the bikes well, mostly it makes the streets and city look more neat.
We are supplying more choices for customers to choose, different countries like different bike racks including double deck bike rack. Take below Japanese style double decker bike rack, Japanese like it more than other countries, you can check following pictures of our new design double decker bike rack if it is your taste.
It is very easy to pull down, so it turns out to be easy for elderly, girls, besides the lower racks can slide in both directions.