# Improved Deep Spectral Convolution Network For Hyperspectral Unmixing With
Multinomial Mixture Kernel and Endmember Uncertainty

Savas Ozkan, and Gozde Bozdagi Akar S. Ozkan and G.B. Akar are with the
Department of Electrical and Electronics Engineering, Middle East Technical
University, 06800, Ankara, Turkey.

###### Abstract

In this study, we propose a novel framework for hyperspectral unmixing by
using an improved deep spectral convolution network (DSCN++) combined with
endmember uncertainty. DSCN++ is used to compute high-level representations
which are further modeled with Multinomial Mixture Model to estimate abundance
maps. In the reconstruction step, a new trainable uncertainty term based on a
nonlinear neural network model is introduced to provide robustness to
endmember uncertainty. For the optimization of the coefficients of the
multinomial model and the uncertainty term, Wasserstein Generative Adversarial
Network (WGAN) is exploited to improve stability and to capture uncertainty.
Experiments are performed on both real and synthetic datasets. The results
validate that the proposed method obtains state-of-the-art hyperspectral
unmixing performance particularly on the real datasets compared to the
baseline techniques.

###### Index Terms:

Hyperspectral Unmixing, Deep Spectral Convolution Network, Endmember
Uncertainty, Spectral Variability, Multinomial Mixture Model, Wasserstein
Generative Adversarial Model

##  I Introduction

A hyperspectral imagery can provide comprehensive knowledge about the Earth
surface in the form of narrow-band spectra that has been used in the variety
of remote sensing applications. Indeed, the separation of constituent
materials (i.e., endmembers) and their fractions (i.e., abundances) which
contribute to the measured hyperspectral data is a fundamental step since the
spectra of the materials can be subjected to complex interactions [1, 2]. As a
result, they can be mixed in different fractions due to the limitations of
hyperspectral sensors and aggravate the problem to be overcome under these
scattering effects. Spectral unmixing aims to determine the fractions of the
materials blindly from the mixture data that enables to analyze the data more
easily [1, 3, 4, 5].

Even if the mixture data can be simplified with a linear formulation (i.e.,
linear mixing model) [6], in real applications, the model remains weak and
underestimates the solution, in particular for the cases such as atmospheric /
lighting conditions and variability of the intrinsic material spectra. To
straighten the model, nonlinear models reformulate the solution by employing
nonlinearity and sparsity in their assumptions [6, 7, 8, 9, 10, 11]. In the
literature, the bilinear model is a commonly used technique for handling
nonlinear scenarios where the bilinear interactions of linear abundances are
taken into account in the model. As stated [6], the methods based on the
bilinear model mainly differ with the additional parameter terms and
constraints exploited in their formulations.

However, both models assume that a hyperspectral scene is composed of a
limited number of fixed endmember spectra. Hence, the performance of the
methods for endmembers and their abundances may be limited because the
possible spectral variability / endmember uncertainty resulting from the data
is not taken into account. Fig. 1 illustrates the spectra of pixels from the
Jasper and Urban datasets [12, 13] that are labeled as pure materials. It can
be seen that each material can form differently even for the same scene under
the assumption of spectral variability. For this purpose, the conventional
methods suffer from the disadvantages of the models (i.e., linear and non-
linear models), and the representations are degraded by the use of false end
members and their combinations. More recently, the angler-based loss function
methods have achieved state-of-the-art performance that can partially overcome
their limitations [2]. Similarly, the idea of stack autoencoders is adapted to
learn blindly spectral signatures. This provides robustness for outliers and
noise sensitive data due to the large set of trainable parameters [14].
Moreover, the methods [2, 15] explain the importance of sparsity for
hyperspectral data. [16] proposes a deep convolution network to estimate
abundances while using both spectral and spatial neighborhoods. However, due
to the lack of examples, the solution is overfitted.

In addition to the linear and nonlinear models, a variety of techniques
consider the spectral variability in their formulations. In general, the
definition of spectral variability assumption is simply a linear/non-linear
combination of true endmembers $\mathbf{e}_{k}\in\rm I\\!R^{D}$ and abudances
$y_{k}\in\rm I\\!R^{1},k=1,...,K$ for each pixel $\mathbf{x}\in\rm I\\!R^{D}$
so that there is a random uncertainty term $\mathbf{\gamma}_{k}\in\rm
I\\!R^{D}$ that also contributes to the material response [17]:

$\displaystyle\mathbf{x}=\sum_{k=1}^{K}{(\mathbf{e}_{k}+\mathbf{\gamma}_{k})y_{k}}+\eta,\hskip
5.69054pts.t.\hskip 3.41432pty_{k}\geq 0,\hskip
3.41432pt\sum_{k=1}^{K}{y_{k}}=1$ (1)

where $K$ is the number of endmembers exhibited in the scene and $D$ is the
spectral bands of a pixel. Also, $\eta\sim\mathcal{N}(0,1)$ is the additive
Gaussian noise to account the possible noise sources. Note that
$\mathbf{\gamma}_{k}$ can vary for each pixel and material, regardless of the
scene, which virtually defines the uncertainty of the data.

In the literature, the studies that intend to be robust to the endmember
uncertainty are based on two different sets of assumptions. One of them
exploits a library of on-hand material spectra (i.e., the set-based
assumption) so as to determine abundances [18, 19, 20, 21, 22] while the other
uses the statistical distributions of the materials (i.e., the distribution-
based assumption) [23, 24, 25, 17]. For the set-based assumption, the
straightforward technique to represent the spectral variability is the
multiple endmember spectral mixture analysis (MESMA) [18]. This algorithm
intuitively unmixes data iteratively by trying all possible material spectra
presented in the library until the error/conditions are completely satisfied
or reached to the desired values. The main drawback is that this method
requires a lot of computational effort, especially when a large body of a
spectra library is used. It is inevitable to conduct a brute-force search for
every possible models is intractable. Its variants have an objective to reduce
the computation complexity and makes the assumption more tractable. Most
variants in the literature focus on reducing library size by selecting the
representative spectra from the library [20, 21, 22]. The rest [19] modifies
the search algorithm which incrementally determines the possible solutions
from a subset of endmembers so that the complexity exhausted in the search
step is moderately reduced.

Similarly, the distribution-based models present the final endmembers as if
they are from some distributions, so that the flexibility for the spectral
variability can be improved. For this purpose, different distribution types
can be used in the literature to represent the endmembers (i.e., normal
composition model (NCM) [23] or beta composition model (BCM) [25]). However,
Zhou et. al. [17] states that the true representation of endmembers might not
be approximated to an unimodal distribution and the mixture of sub-
distributions (i.e., Gaussian Mixture Model (GMM)) may be more appropriate for
modeling real data in particular. It should be noted, however, that the method
needs high memory requirements which makes the algorithm impractical for any
real scene, even at medium spatial resolutions. Furthermore, even if the
expectation maximization is one of the prominent techniques for optimizing GMM
coefficients [26], it suffers from convergence to global solutions and it is
not reliable to learn the ideal multinomial distributions as explained in [26,
27, 28, 29, 30].

In this study, we extend the Deep Spectral Convolution network (DSCN)
framework which is combined with the endmember uncertainty and the multinomial
mixture kernel for hyperspectral unmixing. For this purpose, improved DSCN
(i.e., DSCN++) is exploited to compute high-level representations. Later,
these representation are modeled with Multinomial Mixture Model to estimate
abundance maps. In addition, a new trainable uncertainty term based on a
nonlinear neural network model is introduced in the reconstruction step to
provide robustness to endmember uncertainty. For the optimization of the
parameters of the multinomial model and the uncertainty terms, Wasserstein
Generative Adversarial Network (WGAN) is used to improve stability and to
capture the uncertainty. Note that use of Wasserstein Generative Adversarial
Network (WGAN) for the first time is fairly novel compared to baseline
techniques due to its advantages and performance. To this end, hyperspectral
data is unmixed and the abundances are estimated blindly from a set of
precomputed endmember spectra in an end-to-end learning scheme.

In summary, the overall model takes hyperspectral data as input and computes
low-dimensional representations by using Improved DSCN (DSCN++). Then, these
representations are modeled with the Multinomial Mixture Model to estimate
abundance maps. In the reconstruction step, abundance maps are fed to two NN
modules along with the projection of endmember signatures. To this end,
individual outputs are summed and the input is reconstructed.

The contributions presented in the manuscript can be summarized as follows:

  * •

First, we improve the Deep Spectral Convolution Network (DSCN) model [31] by
introducing modifications to the architecture. This practically reduces the
chance of vanishing gradient problem and enhances the selectivity of the
filters. To this end, the trainable parameters tend to converge to a more
stable solution which is one of the drawbacks highlighted in the baseline
algorithm [31].

  * •

Second, we introduce a multinomial mixture kernel based on a neural network
(NN) architecture that estimates the abundances per-pixel by mimicking the
mixture of Gaussian distributions. To improve the effectiveness, the latent
features obtained from the DSCN are used as the input for the multinomial
mixture kernel. In particular, the parameter set used for the kernel is
optimized with the Wasserstein GAN [32, 33] for the first time which leads to
more accurate and stable solutions than the expectation maximization, the
Kullback-Leibler (KL) [34] divergence or the mean square/absolute error (ME)
[35]. Moreover, WGAN model mitigates the limitations of EM-like methods such
as sensitivity for parameter initialization and discontinuity that reduces the
performance of GMM models. To this end, a flexibility of the model is defined
with robustness to spectral variability.

Figure 1: Spectra of pixels from the Jasper (first row) and Urban (second row)
datasets that correspondence to two different pure materials. It can be seen
that hyperspectral data can be highly influenced by the assumption of spectral
variability.

  * •

Third, we define a new trainable uncertainty term based on a nonlinear NN
model which simply provides robustness to the assumption of endmember
uncertainty for the model. To do so, the estimated abundances are conditioned
with an additional noise distribution as in [32, 33] to synthesize possible
uncertainty. Furthermore, the similar architecture is exploited to minimize
the residue estimated from precomputed endmember spectra in the reconstruction
step. To this end, this auxiliary model practically helps to converge to a
better solution for abundance estimation.

  * •

Note that all these modifications are formulated as a full hyperspectral
unmixing pipeline with NN modules and it is optimized by a stochastic
gradient-based solver in an end-to-end learning scheme.

The rest of the manuscript is structured as follows. First, we define the
concept of the DSCN and explain the modifications presented in the manuscript.
Then, the details of the multinomial mixture kernel and endmember uncertainty
are described. Later, the optimization step for learning trainable parameters
are summarized. Lastly, the experimental results obtained on several real and
synthetic datasets are presented and we conclude the manuscript.

Figure 2: Flow of the proposed method. This model takes hyperspectral data
$\mathbf{x}$ as input and computes high-level representations $\mathbf{z}$ by
using Improved DSCN. Then, these representations are modeled with Multinomial
Mixture Model (i.e., optimized by WGAN) to estimate abundance maps
$\mathbf{\hat{y}}$. In the reconstruction step, abundance maps are fed to two
NN modules with the multiplication of endmember signatures. To this end,
individual outputs are summed and the input is reconstructed.

##  II Improved Deep Spectral Convolution Network

In this section, we first provide a basic definition related to the Deep
Spectral Convolution Networks (DSCNs). Then, the modifications for the
architecture as well as the multinomial mixture kernel and the endmember
uncertainty are explained in detail. Lastly, the optimization step is
described. To ease the understandability, the flow of the proposed method is
illustrated in Fig. 2.

The original DSCN architecture [31] is ultimately equivalent to the standard
autoencoder formulation where it initially encodes the input spectra to a
hidden latent feature $\mathbf{\hat{y}}=Enc(\mathbf{x};\theta_{e})$ and then
decodes this latent feature to recompute the reconstructed version of the
original spectra $\mathbf{\hat{x}}=Dec(\mathbf{\hat{y}};\theta_{d})$.
Ultimately, iterative transformations are performed and the trainable
parameter sets for both encoder $\theta_{e}$ and decoder $\theta_{d}$ parts
are updated based on the reconstruction error calculated by a loss
function(s). To this end, $\mathbf{\hat{y}}$ corresponds to the abundances
per-pixel.

In particular, the main objective of the deep spectral convolution networks
(DSCNs) is that on the contrary to linear layers [2], it blindly extracts
discriminative latent features over data by abstracting the local spectral
characteristics. This is important because the dimensionality of the data
limits the model to achieve an ideal solution through affine transformations
(i.e., corresponds to matrix operations as explained in [2]) learned from
unsupervised data. This can lead to irregularities / holes in the course of
partitioning the feature space. To overcome this limitation, either supervised
data [36, 37] or splitting the data into several overlapping/non-overlapping
parts [38] is needed to be utilized (the detailed information about this
discussion can be found in [36, 31]). Note that the assumption that the data
is split into multiple overlapping / non-overlapping portions directly shares
similar targets in the DSCNs with an unsupervised setting.

Furthermore, the DSCN architecture is mainly composed of 1D stacked
convolutions (i.e., with normalization, non-linear activation and pooling
layers) per-pixel to reveal the latent spectral characteristics of data.

TABLE I: Architecture of the improved DSCN. The notations of Conv1D(W, K, S) indicates 1D convolution whose output dimension, kernel size and stride are W, K and S respectively. Furthermore, SN denotes Spectral Normalization, BN indicates Batch Normalization and AvgP(K) corresponds to Average Pooling layer with kernel size and stride of K.  Index | Index of Inputs | Operation(s) | Output Shape
---|---|---|---
(1) | - | Input Spectra | D
(2) | (1) | Conv1D(10, 21, 1) | D $\times$ 10
(3) | (2) | LReLU + AvgP(5) + BN | D/5 $\times$ 10
(4) | (3) | Conv1D(10, 3, 1) | D/5 $\times$ 10
(5) | (3) | Conv1D(10, 5, 1) | D/5 $\times$ 10
(6) | (3) | Conv1D(10, 7, 1) | D/5 $\times$ 10
(7) | (4),(5),(6) | Concat | D/5 $\times$ 30
(8) | (7) | LReLU + AvgP(2) + SN | D/10 $\times$ 30
(9) | (8) | Conv1D(10, 3, 1) | D/10 $\times$ 10
(10) | (9) | LReLU + AvgP(2) + SN | D/20 $\times$ 10
(11) | (10) | Linear(M) + LReLU | M

However, the lack of parameter convergence and parameter initialization are
some issues in order to obtain stable results for the DSCNs [31]. Therefore,
the architecture should be restructured so that the error calculated from the
loss function is effectively propagated to the earlier layers in the model,
otherwise the overfitting risk for the architecture can be multiplied. Note
that this problem is renowned as the vanishing gradient problem in the
literature [37, 40]. To this end, the modified architecture for
$\mathbf{z}=H(\mathbf{x};\theta_{e})$ is illustrated in Table I.

For the clarity and understandability of the terminology, we divide the
encoder part $Enc(\mathbf{x};\theta_{e})$ into two consecutive steps. First,
the proposed DSCN module is applied to extract high-level latent features
$\mathbf{z}=H(\mathbf{x};\theta_{e})$ (i.e. DSCN++) where $\mathbf{z}\in\rm
I\\!R^{M}$ is the latent feature for each input spectra $\mathbf{x}$. Here,
$M$ denotes the dimension of the latent feature and it can be tuned by the
users (It is set to $10$ throughout the experiments). Then, these latent
features $\mathbf{z}$ are used to estimate the abundances
$\mathbf{\hat{y}}=G(\mathbf{z};\theta_{e})$ (i.e. Multinomial Mixture Kernel
with WGAN) by exploiting the proposed multinomial mixture kernel for the
second step.

For this purpose, initially, the parameter-free ReLU activation (called ReLU)
is replaced by leaky rectified units (LReLU) [41, 42] in the model. This
function is well studied in the literature, especially for regression
problems, as this unit allows small yet non-zero gradients for negative
responses of an input, even if it is not completely active. Moreover, the
influence/resolution of negative responses is adjusted by a constant value
determined as [41].

Second, the order of normalization layers (i.e. spectral normalization [31] or
batch normalization [43]) and spectral convolutions are reversed to reduce the
chance of getting small/no activation(s). Note that spectral normalization
rescales the input by merely leveraging spectral statistics (i.e., mean and
variance of spectral responses of data than batch characteristics). Therefore,
even if the normalized spectral values improve the selectivity of the
representation particularly for the problem, it has also chance to increase
the sparsity with the combination of ReLU [31] (i.e., the chance of vanishing
gradient problem). In Table II, the percentage of active responses in the
model by using the normalization layers before (ours) and after convolution
operations on different datasets. It is clear that the proposed structure
(i.e., pre-normalization) allows more information to forward to the next
layers.

TABLE II: Percentage of active responses in the model for Jasper and Urban datasets. Methods | Jasper | Urban
---|---|---
Post-Normalization | 42.93 | 39.18
Pre-Normalization (ours) | 51.79 | 48.42

In addition, we observe that after initial convolution layer, applying batch
normalization than spectral normalization leads to more accurate/stable
representations since batch normalization can regularize the network as
intended in Dropout layer [40] that reduces the overfitting of the trainable
parameters.

Figure 3: The data distributions before (first column) and after (second
column) applying the improved DSCN to the input spectra for the Jasper (first
row), Urban (second row) and Synthetic (third row) datasets. As seen, the
latent feature space can be more separable than the original material spectra.

When we consider the first convolution layer of the model as a low-level
feature extractor, the consecutive layers aims to compute higher / distinct
representations about data. Also note that the dimensionality of data various
for different sensors. Inspired from [45], the inception module is adapted to
enrich the capacity of the representation. In this way, multiple filters with
different kernel sizes are performed for the same input and filters can select
the best filters by reweighing the responses according to their importance for
the data.

In addition, average pooling layers are leveraged to reduce the dimensionality
of data at each layer. Note that the purpose of using an average operation
instead of a max operation is to reduce the chance of overfitting to the
constantly responding filters [40].

Note that WaveNet model applied for the synthesis of audio signals [46] has a
shared notion as in the proposed model and proves the success for 1D signals.

To demonstrate the effectiveness of the improved DSCNs, Figure 3 visualizes
the data distributions projected to 2D spaces before and after applying the
improved DSCN. Note taht colors are discriminative based on their abundance
values and the PCA is used to reduce the dimensionality. It can be seen that
the improved DSCN transforms the representation of data such that the latent
features become more representative and the materials can be separable from
the others.

In the following section, we will explain the multinomial mixture kernel with
WGAN loss that calculates the abundances per-pixel by exploiting hidden latent
features obtained by the improved DSCN.

###  II-A Multinomial Mixture Kernel For Abundance Maps

After high-level representations $\mathbf{z}$ are computed from data
$\mathbf{x}$, corresponding abundances $\mathbf{y}$ need to be obtained in an
unsupervised manner. For this purpose, a NN-based multinomial mixture model is
utilized.

For the assumption of endmember uncertainty, compositional models sampled from
an unimodal distribution estimate the distributions as follows:

$\displaystyle
p(y_{k}|\mathbf{z})=\mathcal{N}(\mathbf{z}|\mathbf{\mu}_{k},\,\mathbf{\Sigma}_{k}),$
(2)

Here, abundance maps $\mathbf{y}$ are estimated from given $\mathbf{z}$ by
using a set of Gaussian distributions that are composed of mean
$\mathbf{\mu}_{k}\in\rm I\\!R^{M}$ and covariance matrix
$\mathbf{\Sigma}_{k}\in\rm I\\!R^{M\times M}$ where $k=1,...,K$.

As previously mentioned, these kernels can be misleading since the data cannot
be directly approximated to an unimodal distribution [17] and multinomial
mixture models need to be utilized to model the data properly.

For this purpose, Eq. 2 is generalized with Gaussian Mixture Model (GMM) by
computing the distribution $p(y_{k}|\mathbf{z})$ from the statistical
perspective:

$\displaystyle
p(y_{k}|\mathbf{z})=\sum_{n\in\mathcal{K}}p(n)p(y_{k}|\mathbf{z},n),$ (3)

Here, $n$ is a random variable whose sample space is $\mathcal{K}$. Moreover,
the probabilistic values correspond to $p(n)=\pi_{k,n}$ and
$p(y_{k}|\mathbf{z},n)=\mathcal{N}(\mathbf{z}|\mathbf{\mu}_{k,n},\,\mathbf{\Sigma}_{k,n}),n=1,...,N,k=1,...,K$.
Hence, the model formulates the estimation of abundances such that each
material abundance $y_{k}$ can be represented with the mixture of multivariate
Gaussian distributions for a latent feature $\mathbf{z}$ as follows:

$\displaystyle
p(y_{k}|\mathbf{z})=\sum_{n=1}^{N}\pi_{k,n}\mathcal{N}(\mathbf{z}|\mathbf{\mu}_{k,n},\,\mathbf{\Sigma}_{k,n})\hskip
0.28453pt.$ (4)

where $N$ is the number of components in the GMM and the constraints must be
satisfied $\pi_{k,n}\geq 0,\sum_{n=1}^{N}\pi_{k,n}=1$. Moreover, $N\geq K$
assures that endmembers are overly represented with a mixture of latent
Gaussian components that practically provides robustness for the endmember
uncertainty. Moreover, $\pi_{k,n}$ is the weight term which determines the
influence of the $n^{th}$ GMM component to the $k^{th}$ abundance value.

Inspired by methods based on multinomial distributions, we claim that the
distributions (i.e., components) can be imposed on abundances so that each
abundance for a material can be overly expressed with the distributions. Note
that this is ultimately related to the idea of spectral variability. For this
purpose, novel NN kernel and optimization function are proposed so that it
computes the material abundances from the latent features $\mathbf{z}$ (i.e.
Multinomial Mixture Kernel).

For our model, we mimic the GMM with a neural network (NN) module and a novel
multinomial mixture kernel is proposed. This model is trainable in an end-to-
end manner and necessary parameters in the multinomial mixture kernel can be
learned by error backpropagation. Note that a novel optimization function
based on WGAN will be introduced in the subsection II-B

In particular, the distance between the latent representation $\mathbf{z}$ and
$n^{th}$ distribution in Eq. 4 can be calculated simply by the multivariate
mahalanobis distance:

$\displaystyle
d_{k,n}=(\mathbf{\mu}_{k,n}-\mathbf{z})^{T}\mathbf{\Sigma}_{k,n}^{-1}(\mathbf{\mu}_{k,n}-\mathbf{z}).$
(5)

By expanding the formulation in Eq. 5, it can be roughly written as
$d_{k,n}=\mathbf{z}^{T}\mathbf{\Sigma}^{-1}_{k,n}\mathbf{z}-2\mathbf{\mu}_{k,n}^{T}\mathbf{\Sigma}^{-1}_{k,n}\mathbf{z}+\mathbf{\mu}_{k,n}^{T}\mathbf{\Sigma}^{-1}_{k,n}\mathbf{\mu}_{k,n}$.
Note that the formulation can be simplified with a linear neural network layer
(i.e., with trainable weights $\mathbf{w}^{1}_{k,n}$ and bias
$\mathbf{b}^{1}_{k,n}$) by losing the form of a quadratic function as
$d_{k,n}=\mathbf{w}^{1}_{k,n}\mathbf{z}+\mathbf{b}^{1}_{k,n}$ where
$\mathbf{w}^{1}_{k,n}\mathbf{z}=\mathbf{z}^{T}\mathbf{\Sigma}^{-1}_{k,n}\mathbf{z}-2\mathbf{\mu}_{k,n}^{T}\mathbf{\Sigma}^{-1}_{k,n}\mathbf{z}$
and
$\mathbf{b}^{1}_{k,n}=\mathbf{\mu}_{k,n}^{T}\mathbf{\Sigma}^{-1}_{k,n}\mathbf{\mu}_{k,n}$.
Furthermore, the distances for all endmembers is normalized with a sigmoid
activation to produce a probabilistic relation (Please note that softmax
distribution also works seamlessly/similarly for this model). To this end, the
normal distribution is resembled with a cumulative distribution function as
follows:

$\displaystyle\mathcal{N}(\mathbf{z}|\mathbf{\mu}_{k,n},\,\mathbf{\sigma}_{k,n})=\frac{1}{1+e^{\mathbf{w}^{1}_{k,n}\mathbf{z}+\mathbf{b}^{1}_{k,n}}}.$
(6)

Note that the kernel is implemented by directly using $\mathbf{\mu}_{k,n}$ and
$\mathbf{\Sigma}_{k,n}$ coefficients (i.e. as a quadratic form) rather than
parameterizing them with a linear NN layer. However, we observed that even if
this model yields stable results (i.e., small variations), the overall
accuracy is low. As explained in [47], this issue can be summarized by the
fact that parameterizing the formulation with a linear NN layer makes the
operations differentiable which enables to obtain better performance. However,
the capacity of representing true relations (i.e., especially for the
covariance matrix) is also underestimated. Therefore, parameterized
coefficients are used through the manuscript.

Similarly, the weight term $\pi_{k,n}$ is calculated by a NN model and the
sum-to-one constraint is supplied with a softmax activation:

$\displaystyle\pi_{k,n}=\frac{e^{\mathbf{w}^{2}_{k,n}\mathbf{z}+b^{2}_{k,n}}}{\sum_{n=1}^{N}e^{\mathbf{w}^{2}_{k,n}\mathbf{z}+b^{2}_{k,n}}}.$
(7)

To this end, an abundance per-material is modeled as in Eq. 4 by replacing the
original components with the differentiable NN layers formulated with Eq. 6
and Eq. 7. To assure the sum-to-one constraint for abundances of each pixel,
L1-norm is applied at the end [48].

Last but not least, note that the proposed kernel decreases the drawback of
NNs to noise, since multinomial mixture models can effectively capture the
uncertainty in the data so that the trainable parameters yield more stable
results as stated in [49, 50].

###  II-B Parameter Learning

To minimize the reconstruction error, we exploit the optimization scheme
proposed in [2]. For this purpose, the Spectral Angle Distance
$C(\mathbf{x},\mathbf{\hat{x}})$ between the input spectra $\mathbf{x}$ and
the reconstructed spectra $\mathbf{\hat{x}}$ is maximized with the KL-
divergence. Note that the similar notion (i.e., maximize the cross entropy) is
used in GAN optimization:

$\displaystyle\begin{aligned} \mathcal{L}_{re}=-\mathbb{E}[KL\big{(}\hskip
0.28453pt1.0||C(\mathbf{x},{\mathbf{\hat{x}}})\hskip
0.28453pt\big{)}]+\lambda_{0}\mathbb{E}[||\mathbf{x}-\mathbf{\hat{x}}||_{1}]\\\
+\lambda_{1}\|\mathbf{\hat{y}}\|_{1}+\lambda_{2}\|\theta_{e}\|_{2}\end{aligned}$
(8)

where $\lambda_{1}$ determines the influence of the $l1$ regularization term
(i.e., sparsity) and $\lambda_{2}$ is used for regularizing the parameters of
the encoder module (except $\mathbf{w}^{1}$s, $\mathbf{b}^{1}$s) with $l2$
norm. In this work, $\lambda_{1}$ and $\lambda_{2}$ are set to $0.4$ and
$10^{-5}$ respectively. Furthermore, $\lambda_{0}$ defines the ratio of the
mean square/absolute error (ME) contributed to the reconstruction error. Note
that ME is particularly critical for synthetic datasets since the angler
distinctiveness is negligible compared to real data due to the additive noise
(we will explain this statement in the experiment section in detail). Hence,
in the default setting, $\lambda_{0}$ is set to 0. Lastly, $\mathbb{E}[.]$
denotes the expectation value of the possible values for the mini-batch.

In particular, the uncertainty exhibited from data cannot be captured with SAD
and Euclidean-based loss functions. In the proposed method, it is modeled with
WGAN loss by providing theoretical and mathematical explanation.

First, the expectation maximization (EM) is the most popular technique to
estimate GMM parameters in the literature and it theoretically guarantees to
obtain a stationary solution by minimizing the negative log-likelihood
function. However, the log-likelihood function highly suffers from to converge
to a global solution by random parameter initialization and this assumption is
highlighted in various studies in the literature. As explained in [26, 27, 28,
29, 30], this stems from the fact that solving negative log-likelihood
function is similar to using the Kullback-Leibler (KL)-divergence and it has
weakness that can be summarized as twofold. First, since GMM is continuous and
locally Lipschitz [27], the optimization method should also be continuous and
differentiable, although this cannot be fully accomplished by the KL-
divergence. Second, the KL-divergence suffers from the local minimum due to
the discontinuity of the function which becomes sensitive to the choice of
initial parameters. Further discussions related to the KL-divergence can be
found in [27].

TABLE III: Architecture of Discriminative module. Differently, batch normalization and strides are used to obtain latent representations. Index | Index of Inputs | Operation(s) | Output Shape
---|---|---|---
(1) | - | Input Spectra | D
(2) | (1) | Conv1D(5, 21, 5) | D/5 $\times$ 5
(3) | (2) | BN + PReLU | D/5 $\times$ 5
(4) | (3) | Conv1D(10, 5, 2) | D/10 $\times$ 10
(5) | (4) | BN + PReLU | D/10 $\times$ 10
(6) | (5) | Conv1D(20, 5, 2) | D/20 $\times$ 20
(7) | (6) | BN + PReLU | D/20 $\times$ 20
(8) | (7) | Linear(5) | D/20 $\times$ 5

Similarly, for various studies regardless of their applications, the
shortcomings of ME for the regression problems are also stated [35]. In short,
the ME-based solutions yield overly smoothed results due to the averaging of
all possible solutions in the manifold space so that the resultant outputs can
be underestimated. Possible solutions for an inverse problem can be
exponentially multiplied if only pixel-wise relations are used. Note that this
assumption is built for the assessment of image quality [51, 52] and due to
its weakness, ME cannot preserve the structural similarity of data.

As stated, the loss function utilized in the proposed method is based on the
KL-divergence and the ME by minimizing the reconstruction error. Inevitably,
the drawbacks of the techniques as we previously discussed and uncertainty
exhibited from data adversely affect the performance. For this purpose, we
exploit the Wasserstein GAN [32, 33] to optimize the trainable parameters
especially for the mixture model and the uncertainty terms.

Briefly, generative adversarial networks (GANs) have shown impressive
performance for learning data distributions in various machine learning (ML)
problems by their ability to transform a known distribution $p$ to an unknown
data distribution $q$. However, the original GANs are potentially not
continuous to minimize the error as explained in [32] and the similar
limitations can be observed as if the KL-divergence is used (i.e.,
discontinuity). Therefore, the Wasserstein GAN is exploited in the proposed
model which has several theoretical benefits over the original GAN
architecture [32, 53, 54].

Formally, Wasserstein distance $W(p,q)$ is constructed to measure the distance
between $p$ and $q$ (generating variables $\mathbf{x_{p}}$ and
$\mathbf{x_{q}}$ respectively.) by using the Kantorovich-Rubinstein duality
[55] and minimizes the distance between two distributions:

$\displaystyle\begin{aligned} \begin{subarray}{c}\\\ max\\\ \|f\|_{L}\leq
1\end{subarray}\bigg{(}\mathbb{E}_{\mathbf{x_{p}}}[f(\mathbf{x_{p}})]-\mathbb{E}_{\mathbf{x_{q}}}[f(\mathbf{x_{q}})]\bigg{)},\end{aligned}$
(9)

Note that $f(.)$ is a 1-Lipschitz function. Inspired by the structure of
Variational AutoEncoder (VAE) [kingma2013auto], this duality formulation can
be adapted to our problem:

$\displaystyle\begin{aligned} \begin{subarray}{c}\\\ max\\\ \|f\|_{L}\leq
1\end{subarray}\bigg{(}\mathbb{E}_{\mathbf{x}}[f(\mathbf{x})]-\mathbb{E}_{\mathbf{\hat{x}}}[\sum_{k=1}^{K}\mathbf{e}_{k}\sum_{n=1}^{N}\pi_{k,n}f(\mathbf{\hat{x}}_{k,n})]\bigg{)},\end{aligned}$
(10)

where $\mathbf{e}_{k}$s and $\pi_{k,n}$ are not random variables. Moreover,
$\mathbf{\hat{x}}_{k,n}$ corresponds to the generated sample from $n^{th}$ GMM
component for each material $k$. Thus, this formulation can be simplified as
follows:

$\displaystyle\begin{aligned} \begin{subarray}{c}\\\ max\\\ \|f\|_{L}\leq
1\end{subarray}\bigg{(}\mathbb{E}_{\mathbf{x}}[f(\mathbf{x})]-\sum_{k=1}^{K}\mathbf{e}_{k}\sum_{n=1}^{N}\pi_{k,n}\mathbb{E}_{\mathbf{\hat{x}}_{k,n}}[f(\mathbf{\hat{x}}_{k,n})]\bigg{)},\end{aligned}$
(11)

In particular, the training process is stabilized by a gradient penalty [33]
so that final loss function for Wasserstein Adversarial model is written as:

$\displaystyle\begin{aligned}
\mathcal{L}_{adv}=\mathbb{E}[D(\mathbf{x};\theta_{adv})]-\mathbb{E}[D(\mathbf{\hat{x}};\theta_{adv})]\hskip
2.84526pt\\\
-\lambda_{pq}\mathbb{E}[(\|\triangledown_{\mathbf{\tilde{x}}}D(\mathbf{\tilde{x}};\theta_{adv})\|_{2}-1)^{2}].\end{aligned}$
(12)

where $D(.)$ is the discriminative module and it is composed of a family of
1-Lipschitz functions [27]. Moreover, $\mathbf{\tilde{x}}$ is sampled
uniformly from $l1$ normalized input and reconstructed spectra (This is
critical since the angler-based metric is used in the model). $\lambda_{pq}$
determines the influence of the gradient penalty [33] for stable parameter
learning.

Note that the Wasserstein distance is continuous and differentiable almost
everywhere. Moreover, compared to the KL-divergence, the loss function is much
smoother so that the method converges to the ideal solution freely from the
initialization. Additionally, as discussed in [35], more accurate results can
be obtained, since it provides robustness to the averaging operations in the
ME.

For the discriminative module $D(.)$, we adapt the terminology of the PatchGAN
which is extensively used in the image translation problem [56]. More
precisely, instead of penalizing the overall structure of spectra with the
discriminative module, the spectral patches of spectra are used for the
purpose of capturing local spectral characteristics of real and reconstructed
spectra. The architecture details are presented in Table III.

###  II-C Spectral Variability Assumption

In Eq. 1, the formulation that encapsulates the spectral variability has, in
comparison to the linear / nonlinear models, an additional term
$\mathbf{\gamma}_{k}$ that defines the uncertainty per-pixel subjected to the
abundances. In particular, this term should also be modeled with an additional
parameter set in the model to obtain state-of-the-art performance, as it
captures the possible deformations of spectra due to the spectral variability
that cannot be accounted by the multinomial mixture kernel.

In the formulation, by simply expanding Eq. 1, it is easy to see that there
are two sets of equations where $\mathbf{e}_{k}\hat{y}_{k}$ denotes the
projection of endmember spectra $\mathbf{e}_{k}$ (i.e., precomputed material
spectra) with the estimated abundances $\hat{y}_{k}$ while the other term
corresponds to $\mathbf{\gamma}_{k}\hat{y}_{k}$ that defines the influence of
the uncertainty for each abundance. For this purpose, we introduce a trainable
NN module $U(\mathbf{\hat{y}},\eta;\theta_{u})$ which takes the estimated
abundances $\mathbf{\hat{y}}$ and Gaussian noise $\eta\sim\mathcal{N}(0,1)$ as
input and generates an additive nonlinear spectra. This spectra ultimately
degrades the reconstructed spectra based on the uncertainty assumption. More
specifically, the uncertainty module $U(.)$ consists of two consecutive
nonlinear NN layers that incrementally increases the dimension of spectra
(i.e., $K+L\rightarrow 20\rightarrow D$) by exploiting the abundance values
and Gaussian noise. Here, $L$ is the dimension of the noise. Final responses
are rescaled with TanH activation as well as a trainable coefficient
$\alpha_{u}\in[0,0.1]$ to limit the influence.

In addition, since the performance of hypersepctral unmixing is highly
influenced by precomputed material spectra $\mathbf{\hat{e}}_{k}$ (i.e., it
can be computed from various endmember extraction techniques), the
representation of endmember spectra can be also improved with the model for
better abundance results. This is reasonable because current state-of-the-art
methods for endmember extraction are based on NN structure [2, 16]. Instead of
optimizing the precomputed material spectra by backpropagation, we have
observed that introducing a new trainable term that refines the material
spectra results in better performance.

For this purpose, the true material spectra is decomposed into
$\mathbf{e}_{k}=\mathbf{\hat{e}}_{k}+\mathbf{r}_{k}$ where
$\mathbf{\hat{e}}_{k}$ is the precomputed material spectra and
$\mathbf{r}_{k}$ is the possible residue obtained by the endmember method.
Similarly, this residue term can be modeled by exploiting the similar
architecture as in the uncertainty with no shared parameters as
$\mathbf{r}_{k}=U(\mathbf{\hat{y}};\theta_{r})$ and $\alpha_{r}\in[0,0.05]$.
Hence, Eq. 1 in the manuscript that corresponds to
$\mathbf{\hat{x}}=Dec(\mathbf{\hat{y}};\theta_{d})$ can be rewritten as
follows:

$\displaystyle\mathbf{\hat{x}}=\sum_{k=1}^{K}{\mathbf{e}_{k}\hat{y}_{k}}+\alpha_{r}U(\mathbf{\hat{y}};\theta_{r})+\alpha_{u}U(\mathbf{\hat{y}},\eta;\theta_{u})$
(13)

In particular, parameters that capture the uncertainty of endmembers must be
optimized with a novel scheme to achieve state-of-the-art performance. In the
following section, we will introduce a novel loss function based on WGAN that
is applied to hyperspectral unmixing problem for the first time.

###  II-D Implementation Details

Since the proposed model consists of different modules, five individual
parameter sets should be learned from data for parameter optimization. These
are $\theta_{g}$ (corresponds to $\mathbf{w}_{k,n}^{1}$s,
$\mathbf{b}_{k,n}^{1}$s), $\theta_{e}$, $\theta_{r}$, $\theta_{u}$ and
$\theta_{adv}$. Therefore, loss functions and corresponding coefficient terms
for each set are optimized as follows:

$\mathcal{L}_{\theta_{g}}=0.01\hskip
1.13809pt\mathcal{L}_{re}+0.1\mathcal{L}_{adv},$ (14a)
$\mathcal{L}_{\theta_{e}}=\mathcal{L}_{re},$ (14b)
$\mathcal{L}_{\theta_{r}}=0.001\hskip 1.13809pt\mathcal{L}_{re},$ (14c)
$\mathcal{L}_{\theta_{u}}=0.001\hskip 1.13809pt\mathcal{L}_{adv},$ (14d)
$\mathcal{L}_{\theta_{adv}}=\mathcal{L}_{adv}.$ (14e)

From Eq. 14, it can be seen that the Wasserstein GAN is used primarily to
optimize the parameters of $\theta_{g}$ and $\theta_{u}$ that implicitly
provides robustness to the parameters for spectral variability / uncertainty.
GAN models can be powerful enough to generate a proper noise space based on
the definition of the problem. This is critical since the uncertainty
exhibited from hyperspectral data can be precisely captured by the Wasserstein
distance. Moreover, Wasserstein distance has several theoretical benefits to
obtain a stable and accurate solution as explained throughout the manuscript.

Moreover, to optimize the latent representations learned from data, the DSCN
parameters are updated by exploiting only the reconstruction loss
$\mathcal{L}_{re}$. In addition, the learning rates for multinomial mixing
kernel and auxiliary functions are scaled with the constants to decrease the
learning speed compared to the DSCN module. This is because DSCN parameters
should converge faster to get more reliable abundances as expected.

Lastly, Adam stochastic optimizer [57] is utilized to update the parameters
for each mini-batch by empirically defining the first-order moment term
$\beta_{1}$ as 0.7. Note that the mini-batch size is fixed to 64. Also, we set
the learning rate and the number of training iterations to 0.002 and 10K
respectively for all datasets. The proposed algorithm is implemented on Python
by extensively leveraging Tensorflow framework111Project Page:
https://github.com/savasozkan/dscn . Also, additional results are also
provided.

##  III Experiments

Experiments to compare the performance of the proposed method (DSCN++) with
the baseline techniques are conducted on one synthetic and two real datasets.
For this purpose, Linear Mixture Model (LMM) [7], Generalized Bilinear Model
(GBM) [58], Post-Nonlinear Mixing Model (PPNM) [59], Multilinear Mixing Model
(MLM) [7] and Normal Compositional Model (NCM) [60] are selected for the
hyperspectral unmixing baselines that are extensively used methods in the
literature (codes are available on the web.). Note that the recent method in
[17] needs high memory requirement which is impossible to compute the
algorithm on real datasets especially by maintaining their original spatial
resolutions.

Furthermore, Spatial Compositional Model (SCM) [61], Distance-MaxD (DMaxD) [3]
and Sparse AutoEncoder Network (EndNet) [2] are exploited to estimate
endmember spectra from the scenes.

To evaluate the unmixing performance and compare the estimated abundances with
the ground truth, Root Mean Square Error (RMSE) is utilized:

$\displaystyle
RMSE(\mathbf{y},\mathbf{\hat{y}})=\sqrt{\frac{1}{L}\|\mathbf{y}-\mathbf{\hat{y}}\|^{2}_{2}}.$
(15)

where $\mathbf{y}$ denotes the true abundance per-pixel while
$\mathbf{\hat{y}}$ indicates the estimated ones by the different methods.
Also, $L$ corresponds to the total number of pixels in the scene.

Figure 4: Endmember spectra (first row) and abundance maps (second row) for
the Synthetic dataset.

In following sections, we first provide the quantitative performance for each
dataset that clearly demonstrates the impact of the proposed method by
highlighting each individual step as well as different parameter
configurations in detail. Note that $N$ is the only coefficient needed to be
tuned by the user. Therefore $N$ is determined for each step based on their
overall results in the experiments (As expected, $N$ (the number of mixture
components) can vary for each dataset as well as each endmember method based
on its distinctiveness). Later, the comparisons with the baselines are
explained. For reliability of assessments, tests are repeated 20 times for the
proposed method, thus mean and standard deviation are reported.

To ease the understandability, some abbreviations used in the Tables can be
summarized as follows: DSCN indicates the deep spectral convolution network as
explained in the section II. WGAN indicates the use of Wasserstein GAN in the
parameter optimization as mentioned in the subsection II-B. Lastly, EU stands
for the usage of NN models for the uncertainty term and the auxiliary model
described in the subsection II-C.

TABLE IV: RMSE performance on the synthetic dataset for different configurations. $\times 10^{-2}$ | DMaxD [3] | SCM [61] | EndNet [2]
---|---|---|---
$N=4$ | 8.32 $\pm$1.4 | 4.47 $\pm$0.7 | 9.20 $\pm$1.4
$N=8$ | 8.15 $\pm$1.4 | 4.36 $\pm$0.9 | 9.14 $\pm$1.0
$N=16$ | 8.18 $\pm$1.0 | 3.97 $\pm$0.5 | 9.06 $\pm$1.1
$N=24$ | 7.55 $\pm$1.1 | 4.42 $\pm$0.6 | 9.51 $\pm$1.2
w$\setminus$o DSCN | 10.78 $\pm$1.2 | 10.96 $\pm$1.5 | 13.37 $\pm$1.9
w$\setminus$o EU | 5.43 $\pm$0.3 | 3.81 $\pm$0.9 | 6.43 $\pm$0.4
w$\setminus$o WGAN | 5.23 $\pm$0.4 | 3.22 $\pm$1.2 | 8.88 $\pm$0.3

###  III-A Synthetic Dataset

Figure 5: Spectra of pixels from the Synthetic dataset. It can be seen that
even if random noise is added in the construction of data, spectral
variability cannot be exhibited truly and it can be misleading for learning-
based methods since learning-based methods (generative-based models in
particular) can learn this random noise easily.

The unmixing methods are tested on the synthetic dataset released for the
experiments in [17]. The dataset consists of four constituent materials from
the ASTER spectral library as limestone, basalt, concrete and asphalt. The
spectral range of data is from $0.4\mu$m to $14\mu$m and the spectral
dimension is resampled to 200. Furthermore, the spatial resolution of the
scene is $60\times 60$. For the scene, limestone is set as the background
material, the rest is randomly placed to the scene by maintaining the
Gaussian-shaped distributions. The endmember spectra and abundance maps for
the scene are illustrated in Fig. 4.

Note that even if synthetic datasets are conventional to test the
hyperspectral unmixing performance, it is not quite possible to capture all
adverse conditions observed in the real data (such as spectral variability).
Hence, it can be hard to validate the performance truly under this strong
bias. To degrade the bias, the only practical solution is to add random noise
on the data which cannot reflect the true physical conditions (i.e., spectral
variability) observed in the real data. More precisely, Fig. 1 and Fig. 5 show
the severe difference between real and synthetic datasets for pure materials.
It can be seen that the angler distinctiveness is quite small (i.e., spectra
of pixels are nearly concentrated on a single response) and use of angler-
based metric can be misleading for the dataset. Hence, only for the synthetic
dataset, the KL-divergence with SAD in the reconstruction loss is replaced
(for SCM) or jointly used (for DmaxD, Endnet) with the ME (i.e.,
$\lambda_{0}=1$).

The experimental results for each individual step and different parameter
configuration settings are reported in Table IV. From the results, the
proposed method achieves significantly better performance without EU and WGAN
steps. This is quite reasonable since the data is degraded by random noise to
simulate the spectral variability and the proposed model ultimately learns the
noise from data than the assumption of spectral variability. Hence, this leads
to the corruption for the model. However, note that true impacts of these
steps can be clearly observed for real datasets. Therefore, further
discussions about these steps are reserved for the real datasets.

The RMSE performance for the baselines are shown in Table V. From the results,
the proposed method obtain slightly better performance compared to the
baselines. Note that the proposed method is effective on the synthetic dataset
even if it does not exploit all its features (i.e., EU and WGAN).

Lastly, note that [17] obtains $1.53\times 10^{-2}$ RMSE performance for this
synthetic data. However, due to its computation requirement and differences
for the endmember estimation, it is not fully compatible with our method.

TABLE V: RMSE performance on the synthetic dataset for baselines. The best results are indicated as Bold-Red while the second best results are shown as Bold. $\times 10^{-2}$ | DMaxD [3] | SCM [61] | EndNet [2]
---|---|---|---
LMM [7] | 5.17 | 3.36 | 17.62
GBM [58] | 5.62 | 3.26 | 17.39
PPNM [59] | 7.66 | 4.01 | 9.34
MLM [7] | 6.52 | 3.85 | 6.71
NCM [60] | 4.09 | 3.85 | 11.17
DSCN++ | 5.23 | 3.22 | 6.43

Figure 6: Endmember spectra (first row) and abundance maps (second row) for
the Urban dataset.

###  III-B Urban Dataset

The spatial resolution of the scene is $307\times 307$. The spectral range
covers from 10 nm to 2500 nm and it is composed of 210 spectral bands. Several
spectral bands (1-4, 76, 87, 101-111, 136-153 and 198-210) are removed from
data due to the water-vapor absorption and atmospheric effects, thus the final
spectral band is reduced to 162. Four material types are observed in the
scene: asphalt, grass, tree and roof [13]. The endmember spectra and abundance
maps for the scene are illustrated in Fig. 6.

The performance for the proposed method under different configurations is
illustrated in Table VI. It can be seen that the selection of $N$ larger than
$K$ improves the performance for all endmember techniques. Note that the
optimum $N$ value can vary since the discriminative power of the
representation can be highly influenced by the endmember extraction techniques
regardless of how spectral variability affects the scene. In addition, the
results confirms that use of DSCN in order to reduce the dimension of data and
to obtain compact representations yields better results. Moreover, for the
models without the DSCN, the variations in the performance are significantly
increased due to the fact that high dimensionality can lead to
irregularities/holes for the space as previously described.

TABLE VI: RMSE performance on the Urban dataset for different configurations. $\times 10^{-2}$ | DMaxD [3] | SCM [61] | EndNet [2]
---|---|---|---
$N=4$ | 15.96 $\pm$2.6 | 13.26 $\pm$0.6 | 8.27 $\pm$0.4
$N=8$ | 15.47 $\pm$1.9 | 13.29 $\pm$0.8 | 8.05 $\pm$0.5
$N=16$ | 16.32 $\pm$3.2 | 13.24 $\pm$0.5 | 8.16 $\pm$0.4
$N=24$ | 16.56 $\pm$2.1 | 13.23 $\pm$0.8 | 8.04 $\pm$0.3
w$\setminus$o DSCN | 22.12 $\pm$2.6 | 16.63 $\pm$4.6 | 11.42 $\pm$0.8
w$\setminus$o EU | 34.14 $\pm$0.3 | 15.91 $\pm$0.3 | 8.38 $\pm$0.4
w$\setminus$o WGAN | 31.43 $\pm$1.6 | 15.76 $\pm$1.0 | 8.52 $\pm$0.5
TABLE VII: RMSE performance on the Urban dataset for baselines. The best results are indicated as Bold-Red while the second best results are shown as Bold.  $\times 10^{-2}$ | DMaxD [3] | SCM [61] | EndNet [2]
---|---|---|---
LMM [7] | 31.04 | 28.00 | 27.59
GBM [58] | 31.04 | 27.00 | 27.08
PPNM [59] | 36.54 | 28.17 | 15.06
MLM [7] | 31.43 | 16.34 | 9.48
NCM [60] | 31.04 | 28.00 | 27.51
DSCN++ | 15.47 | 13.23 | 8.04

Similarly, considering endmember uncertainty with additional NN models
practically improves the performance. For the Wasserstein GAN, it helps to
converge to the ideal solution more likely. Moreover, it achieves more stable
results by which the influence of parameter initialization is decreased in the
model.

The performance is compared with the baselines in Table VII and the proposed
method outperforms all hyperspectral unmixing techniques significantly. In
particular, the usage of WGAN and EU steps improve the performance of DmaxD
such that it leads to a highly accurate solution where the baselines stuck to
stationary results.

Figure 7: Endmember spectra (first row) and abundance maps (second row) for
the Jasper dataset.

###  III-C Jasper Dataset

The Jasper dataset is captured by the AVIRIS sensor and the dimension is
$100\times 100$ pixels. Several channels (1-3, 108-112, 154-166 and 220-224)
are discarded due to the atmospheric effects and water-vapor absorption. The
spectral band is of 198. There are four constituent materials in the scene as
tree, water, soil and road [12]. Fig. 6 shows the endmember spectra and
abundance maps for the scene.

The detailed performance for different configuration settings and steps are
reported in Table VIII. It can be seen that similar results can be observed as
stated for the Urban dataset. A single outlier is that SCM without the DSCN
achieves slightly better performance compared to the baseline. This stems from
the fact that the SCM representation space can be more discriminative than the
DSCN for this dataset.

TABLE VIII: RMSE performance on the Jasper dataset for different configurations. $\times 10^{-2}$ | DMaxD [3] | SCM [61] | EndNet [2]
---|---|---|---
$N=4$ | 12.03 $\pm$0.7 | 18.10 $\pm$0.9 | 6.38 $\pm$0.7
$N=8$ | 12.29 $\pm$0.9 | 18.33 $\pm$1.1 | 6.61 $\pm$0.7
$N=16$ | 12.12 $\pm$0.7 | 17.92 $\pm$0.8 | 6.31 $\pm$0.7
$N=24$ | 12.27 $\pm$0.6 | 18.03 $\pm$1.1 | 6.62 $\pm$1.1
w$\setminus$o DSCN | 12.47 $\pm$0.6 | 17.49 $\pm$1.2 | 6.87 $\pm$2.2
w$\setminus$o EU | 14.78 $\pm$0.4 | 20.26 $\pm$0.4 | 8.36 $\pm$0.5
w$\setminus$o WGAN | 15.63 $\pm$0.3 | 19.03 $\pm$0.4 | 10.15 $\pm$0.5
TABLE IX: RMSE performance on the Jasper dataset for baselines. The best results are indicated as Bold-Red while the second best results are shown as Bold. $\times 10^{-2}$ | DMaxD [3] | SCM [61] | EndNet [2]
---|---|---|---
LMM [7] | 15.31 | 19.78 | 22.60
GBM [58] | 15.98 | 19.95 | 21.91
PPNM [59] | 12.99 | 19.78 | 14.72
MLM [7] | 16.59 | 20.73 | 8.13
NCM [60] | 15.29 | 19.56 | 22.42
DSCN++ | 12.03 | 17.92 | 6.31

For the baselines comparisons, the proposed method obtains state-of-the-art
performance as demonstrated in Tab IX. Compared to [31] which corresponds to
the baseline DSCN, the improved model significantly improves the performance
while small variations in the performance can be observed.

##  IV Discussion and Conclusion

In this manuscript, we introduce a novel neural network (NN) framework for
hyperspectral unmixing. It is mainly composed of deep spectral convolution
networks (DSCNs) by addressing the assumption of spectral variability that is
one of the severe physical conditions frequently observed in real data. More
precisely, in real applications, pixel responses for each material can be
varied due to the atmospheric/illumination conditions instead of fixed
material spectra. Therefore, these conditions should be addressed by an
effective model in order to obtain state-of-the-art unmixing performance.

For this purpose, we present critical contributions throughout the manuscript
as follows: First, we make several modifications in the architecture of the
DSCNs to obtain stationary yet more accurate results. Second, a novel NN layer
is proposed to estimate the abundances per-pixel by leveraging the multinomial
mixture kernel. This kernel particularly provides robustness for the
assumption of spectral variability since the distribution-based assumptions
can be more convenient to handle to these adverse conditions as explained in
the literature. Moreover, an extra NN model is utilized to capture the
uncertainty term in the formulation of the spectral variability. Third, the
Wasserstein GAN is exploited in the parameter optimization which leads to more
accurate and stable performance than the KL-divergence and the ME. Lastly, the
model is formulated as a full pipeline based on a NN method that can be
learned by backpropagation and a stochastic gradient solver.

The experimental results validate that the proposed method outperforms well-
known hyperspectral unmixing baselines by a large margin. In particular, the
methods obtains state-of-the-art performances on the real datasets where
severe spectral variability can be observed.

Complexity: As analyzed in [2, 62], a NN model with the stochastic gradient
algorithm significantly decreases the complexity and memory requirement needed
by the model. More precisely, these parameters are independent from the data
size due to the batch-based learning. Hence, the method is more applicable for
large-scale data.

## Acknowledgment

The authors gratefully acknowledge the support of NVIDIA Corporation with the
donation of GPUs used for this research.

## References

  * [1] Z. Yang, G. Zhou, S. Xie, S. Ding, J.-M. Yang, and J. Zhang, “Blind spectral unmixing based on sparse nonnegative matrix factorization,” _IEEE Transactions on Image Processing_ , vol. 20, no. 4, pp. 1112–1125, 2011.
  * [2] S. Ozkan, B. Kaya, and G. B. Akar, “Endnet: Sparse autoencoder network for endmember extraction and hyperspectral unmixing,” _arXiv preprint arXiv:1708.01894_ , 2017.
  * [3] R. Heylen, D. Burazerovic, and P. Scheunders, “Non-linear spectral unmixing by geodesic simplex volume maximization,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 5, no. 3, pp. 534–542, 2011.
  * [4] J. J. Settle, “On the relationship between spectral unmixing and subspace projection,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 34, no. 4, pp. 1045–1046, 1996.
  * [5] N. Keshava and J. F. Mustard, “Spectral unmixing,” _IEEE signal processing magazine_ , vol. 19, no. 1, pp. 44–57, 2002.
  * [6] N. Dobigeon, Y. Altmann, N. Brun, and S. Moussaoui, “Linear and nonlinear unmixing in hyperspectral imaging,” in _Data Handling in Science and Technology_. Elsevier, 2016, vol. 30, pp. 185–224.
  * [7] R. Heylen and P. Scheunders, “A multilinear mixing model for nonlinear spectral unmixing,” _IEEE transactions on geoscience and remote sensing_ , vol. 54, no. 1, pp. 240–251, 2016.
  * [8] B. Somers, K. Cools, S. Delalieux, J. Stuckens, D. Van der Zande, W. W. Verstraeten, and P. Coppin, “Nonlinear hyperspectral mixture analysis for tree cover estimates in orchards,” _Remote Sensing of Environment_ , vol. 113, no. 6, pp. 1183–1193, 2009.
  * [9] Y. Altmann, N. Dobigeon, and J.-Y. Tourneret, “Bilinear models for nonlinear unmixing of hyperspectral images,” in _Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), 2011 3rd Workshop on_. IEEE, 2011, pp. 1–4.
  * [10] J. M. Nascimento and J. M. Bioucas-Dias, “Nonlinear mixture model for hyperspectral unmixing,” in _Image and Signal Processing for Remote Sensing XV_ , vol. 7477. International Society for Optics and Photonics, 2009, p. 74770I.
  * [11] N. Raksuntorn and Q. Du, “Nonlinear spectral mixture analysis for hyperspectral imagery in an unknown environment,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 7, no. 4, pp. 836–840, 2010.
  * [12] F. Zhu, Y. Wang, S. Xiang, B. Fan, and C. Pan, “Structured sparse method for hyperspectral unmixing,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 88, pp. 101–118, 2014.
  * [13] F. Zhu, Y. Wang, B. Fan, S. Xiang, G. Meng, and C. Pan, “Spectral unmixing via data-guided sparsity,” _IEEE Transactions on Image Processing_ , vol. 23, no. 12, pp. 5412–5427, 2014.
  * [14] Y. Su, J. Li, A. Plaza, A. Marinoni, P. Gamba, and S. Chakravortty, “Daen: Deep autoencoder networks for hyperspectral unmixing,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2019.
  * [15] X.-R. Feng, H.-C. Li, J. Li, Q. Du, A. Plaza, and W. J. Emery, “Hyperspectral unmixing using sparsity-constrained deep nonnegative matrix factorization with total variation,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2018.
  * [16] X. Zhang, Y. Sun, J. Zhang, P. Wu, and L. Jiao, “Hyperspectral unmixing via deep convolutional neural networks,” _IEEE Geoscience and Remote Sensing Letters_ , no. 99, pp. 1–5, 2018.
  * [17] Y. Zhou, A. Rangarajan, and P. D. Gader, “A gaussian mixture model representation of endmember variability in hyperspectral unmixing,” _IEEE Transactions on Image Processing_ , vol. 27, no. 5, pp. 2242–2256, 2018\.
  * [18] D. A. Roberts, M. Gardner, R. Church, S. Ustin, G. Scheer, and R. Green, “Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models,” _Remote Sensing of Environment_ , vol. 65, no. 3, pp. 267–279, 1998.
  * [19] R. Heylen, A. Zare, P. Gader, and P. Scheunders, “Hyperspectral unmixing with endmember variability via alternating angle minimization,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 54, no. 8, pp. 4983–4993, 2016.
  * [20] C. Quintano, A. Fernández-Manso, and D. A. Roberts, “Multiple endmember spectral mixture analysis (mesma) to map burn severity levels from landsat images in mediterranean countries,” _Remote Sensing of Environment_ , vol. 136, pp. 76–88, 2013.
  * [21] D. A. Roberts, P. E. Dennison, M. E. Gardner, Y. Hetzel, S. L. Ustin, and C. T. Lee, “Evaluation of the potential of hyperion for fire danger assessment by comparison to the airborne visible/infrared imaging spectrometer,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 41, no. 6, pp. 1297–1310, 2003.
  * [22] P. Dennison and D. Roberts, “Endmember selection for mapping chaparral species and fraction using multiple endmember spectral mixture analysis,” _Remote Sensing of Environment_ , vol. 87, no. 2-3, pp. 295–309, 2003.
  * [23] B. Zhang, L. Zhuang, L. Gao, W. Luo, Q. Ran, and Q. Du, “Pso-em: A hyperspectral unmixing algorithm based on normal compositional model,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 52, no. 12, pp. 7782–7792, 2014.
  * [24] O. Eches, N. Dobigeon, and J.-Y. Tourneret, “Estimating the number of endmembers in hyperspectral images using the normal compositional model and a hierarchical bayesian algorithm,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 4, no. 3, pp. 582–591, 2010.
  * [25] X. Du, A. Zare, P. Gader, and D. Dranishnikov, “Spatial and spectral unmixing using the beta compositional model,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 7, no. 6, pp. 1994–2003, 2014.
  * [26] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the em algorithm,” _Journal of the royal statistical society. Series B (methodological)_ , pp. 1–38, 1977.
  * [27] S. Kolouri, G. K. Rohde, and H. Hoffman, “Sliced wasserstein distance for learning gaussian mixture models,” _arXiv preprint arXiv:1711.05376_ , 2017\.
  * [28] B. Jian and B. C. Vemuri, “Robust point set registration using gaussian mixture models,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 33, no. 8, pp. 1633–1645, 2011.
  * [29] C. Améndola, M. Drton, and B. Sturmfels, “Maximum likelihood estimates for gaussian mixtures are transcendental,” in _International Conference on Mathematical Aspects of Computer and Information Sciences_. Springer, 2015, pp. 579–590.
  * [30] C. Jin, Y. Zhang, S. Balakrishnan, M. J. Wainwright, and M. I. Jordan, “Local maxima in the likelihood of gaussian mixture models: Structural results and algorithmic consequences,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 4116–4124.
  * [31] S. Ozkan and G. B. Akar, “Deep spectral convolution network for hyperspectral unmixing,” _arXiv preprint arXiv:1806.08562_ , 2018.
  * [32] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein generative adversarial networks,” in _International Conference on Machine Learning_ , 2017, pp. 214–223.
  * [33] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. C. Courville, “Improved training of wasserstein gans,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 5769–5779.
  * [34] J. R. Hershey and P. A. Olsen, “Approximating the kullback leibler divergence between gaussian mixture models,” in _Acoustics, Speech and Signal Processing, 2007. ICASSP 2007. IEEE International Conference on_ , vol. 4. IEEE, 2007, pp. IV–317.
  * [35] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang _et al._ , “Photo-realistic single image super-resolution using a generative adversarial network,” _arXiv preprint_ , 2017.
  * [36] M. Norouzi and D. M. Blei, “Minimal loss hashing for compact binary codes,” in _Proceedings of the 28th international conference on machine learning (ICML-11)_. Citeseer, 2011, pp. 353–360.
  * [37] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [38] H. Jegou, M. Douze, and C. Schmid, “Product quantization for nearest neighbor search,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 33, no. 1, pp. 117–128, 2011.
  * [39] A. Nguyen, J. Yosinski, and J. Clune, “Deep neural networks are easily fooled: High confidence predictions for unrecognizable images,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2015, pp. 427–436.
  * [40] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: A simple way to prevent neural networks from overfitting,” _The Journal of Machine Learning Research_ , vol. 15, no. 1, pp. 1929–1958, 2014.
  * [41] A. L. Maas, A. Y. Hannun, and A. Y. Ng, “Rectifier nonlinearities improve neural network acoustic models,” in _Proc. icml_ , vol. 30, no. 1, 2013, p. 3.
  * [42] K. He, X. Zhang, S. Ren, and J. Sun, “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 1026–1034.
  * [43] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” _arXiv preprint arXiv:1502.03167_ , 2015.
  * [44] T. Salimans and D. P. Kingma, “Weight normalization: A simple reparameterization to accelerate training of deep neural networks,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 901–909.
  * [45] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the inception architecture for computer vision,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2016, pp. 2818–2826.
  * [46] A. v. d. Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, “Wavenet: A generative model for raw audio,” _arXiv preprint arXiv:1609.03499_ , 2016.
  * [47] R. Arandjelovic, P. Gronat, A. Torii, T. Pajdla, and J. Sivic, “Netvlad: Cnn architecture for weakly supervised place recognition,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2016, pp. 5297–5307.
  * [48] B. Hapke, “Bidirectional reflectance spectroscopy: 1. theory,” _Journal of Geophysical Research: Solid Earth_ , pp. 3039–3054, 1981.
  * [49] P. Melchior and A. D. Goulding, “Filling the gaps: Gaussian mixture models from noisy, truncated or incomplete samples,” _arXiv preprint arXiv:1611.05806_ , 2016.
  * [50] A. Ozerov, M. Lagrange, and E. Vincent, “Uncertainty-based learning of gaussian mixture models from noisy data,” 2012.
  * [51] Z. Wang, E. P. Simoncelli, and A. C. Bovik, “Multiscale structural similarity for image quality assessment,” in _The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003_, 2003.
  * [52] Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli _et al._ , “Image quality assessment: from error visibility to structural similarity,” _IEEE transactions on image processing_ , pp. 600–612, 2004.
  * [53] C. Frogner, C. Zhang, H. Mobahi, M. Araya, and T. A. Poggio, “Learning with a wasserstein loss,” in _Advances in Neural Information Processing Systems_ , 2015, pp. 2053–2061.
  * [54] G. Peyré, J. Fadili, and J. Rabin, “Wasserstein active contours,” in _Image Processing (ICIP), 2012 19th IEEE International Conference on_. IEEE, 2012, pp. 2541–2544.
  * [55] C. Villani, _Optimal transport: old and new_. Springer Science & Business Media, 2008, vol. 338.
  * [56] P. Isola, J.-Y. Zhu, T. Zhou, and A. A. Efros, “Image-to-image translation with conditional adversarial networks,” _arXiv preprint_ , 2017.
  * [57] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [58] A. Halimi, Y. Altmann, N. Dobigeon, and J.-Y. Tourneret, “Nonlinear unmixing of hyperspectral images using a generalized bilinear model,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 49, no. 11, pp. 4153–4162, 2011.
  * [59] Y. Altmann, A. Halimi, N. Dobigeon, and J.-Y. Tourneret, “Supervised nonlinear spectral unmixing using a postnonlinear mixing model for hyperspectral imagery,” _IEEE Transactions on Image Processing_ , vol. 21, no. 6, pp. 3017–3025, 2012.
  * [60] D. Stein, “Application of the normal compositional model to the analysis of hyperspectral imagery,” in _Advances in techniques for analysis of remotely sensed data, 2003 IEEE Workshop on_. IEEE, 2003, pp. 44–51.
  * [61] Y. Zhou, A. Rangarajan, and P. D. Gader, “A spatial compositional model for linear unmixing and endmember uncertainty estimation,” _IEEE Transactions on Image Processing_ , vol. 25, no. 12, pp. 5987–6002, 2016.
  * [62] L. Bottou, “Large-scale machine learning with stochastic gradient descent,” in _Proceedings of COMPSTAT’2010_. Springer, 2010, pp. 177–186.


#  $f(T)$ cosmology via Noether symmetry

K. Atazadeh1 and F. Darabi1,2
1Department of Physics, Azarbaijan University of Tarbiat Moallem, Tabriz
53741-161, Iran
2Research Institute for Astronomy and Astrophysics of Maragha (RIAAM), Maragha
55134-441, Iran email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider Noether symmetry approach to find out exact cosmological solutions
in $f(T)$-gravity. Instead of taking into account phenomenological models, we
apply the Noether symmetry to the $f(T)$ gravity. As a result, the presence of
such symmetries selects viable models and allow to solve the equations of
motion. We show that the generated $f(T)$ leads to a power law expansion for
the cosmological scale factor.

##  1 Introduction

The expansion of the universe is currently undergoing a period of acceleration
which is directly measured from observations such as Type Ia Supernovae [1],
[2] cosmic microwave background (CMB) radiation [3, 4], large scale structure
[5], baryon acoustic oscillations [6], and weak lensing [7]. There are two
remarkable approaches to explain the late time acceleration of the universe:
One is to introduce some unknown matters called “dark energy” in the framework
of general relativity (for a review on dark energy, see, e.g., [8, 9]). The
other is to modify the gravitational theory, e.g., $f(R)$ gravity [10, 11].

To extend gravity beyond general relativity, “teleparallelism" could be
considered by using the Weitzenböck connection, which has no curvature but
torsion, rather than the curvature defined by the Levi-Civita connection [12],
[13]. This approach was also taken by Einstein [14]. The teleparallel
Lagrangian density described by the torsion scalar $T$ has been promoted to a
function of $T$, $i.e.$, $f(T)$, in order to account for the late time cosmic
acceleration [15] as well as inflation [16]. This concept is similar to the
idea of $f(R)$ gravity. Various aspects of $f(T)$ gravity have been examined
in the literature [17, 18, 19, 20, 21]. In particular, the presence of extra
degrees of freedom and the violation of local Lorentz invariance as well as
the existence of non-trivial frames for $f(T)$ gravity have been noted [21].
Evidently, more studies on $f(T)$ gravity are needed to see if the theory is a
viable one. For a comprehensive review of the teleparallel gravity, the reader
is referred to [22].

In this paper, we consider a flat FRW space-time in the framework of the
metric formalism of $f(T)$ gravity. Following [23], we set up an effective
Lagrangian in which the scale factor $a$ and torsion scalar $T$ play the role
of independent dynamical variables. This Lagrangian is constructed in such a
way that its variation with respect to $a$ and $T$ yields the correct
equations of motion as that of an action with a generic $f(T)$ mentioned
above. The form of the function $f(T)$ appearing in the modified action is
then found by demanding that the Lagrangian admits the desired Noether
symmetry [24, 25]. For a study of the Noether symmetry in $f(R)$ cosmology see
[26]. Similarly, we shall see that by demanding the Noether symmetry as a
feature of the Lagrangian of the cosmological model under consideration, we
can obtain the explicit form of the function $f(T)$. Since the existence of a
symmetry results in a constant of motion, we can integrate the field equations
which would then lead to a power law expansion of the universe.

##  2  $f(T)$ gravity and cosmology

To consider teleparallelism, one employs the orthonormal tetrad components
$e_{A}(x^{\mu})$, where an index $A$ runs over $0,1,2,3$ to the tangent space
at each point $x^{\mu}$ of the manifold. Their relation to the metric
$g_{\mu\nu}$ is given by

$g_{\mu\nu}=\eta_{AB}e^{A}_{\mu}e^{B}_{\nu}\,,$ (1)

where $\mu$ and $\nu$ are coordinate indices on the manifold and also run over
$0,1,2,3$, and $e_{A}^{\mu}$ forms the tangent vector on the tangent space
over which the metric $\eta_{AB}$ is defined.

Instead of using the torsionless Levi-Civita connection in General Relativity,
we use the curvatureless Weitzenböck connection in Teleparallelism [12], whose
non-null torsion $T^{\rho}_{\verb| |\mu\nu}$ and contorsion
$K^{\mu\nu}_{\verb| |\rho}$ are defined by

$\displaystyle T^{\rho}_{\verb| |\mu\nu}\equiv
e^{\rho}_{A}\left(\partial_{\mu}e^{A}_{\nu}-\partial_{\nu}e^{A}_{\mu}\right)\,,$
(2) $\displaystyle K^{\mu\nu}_{\verb|
|\rho}\equiv-\frac{1}{2}\left(T^{\mu\nu}_{\verb| |\rho}-T^{\nu\mu}_{\verb|
|\rho}-T_{\rho}^{\verb| |\mu\nu}\right)\,,$ (3)

respectively. Moreover, instead of the Ricci scalar $R$ for the Lagrangian
density in general relativity, the teleparallel Lagrangian density is
described by the torsion scalar $T$ as follows

$T\equiv S_{\rho}^{\verb| |\mu\nu}T^{\rho}_{\verb| |\mu\nu}\,,$ (4)

where

$S_{\rho}^{\verb| |\mu\nu}\equiv\frac{1}{2}\left(K^{\mu\nu}_{\verb|
|\rho}+\delta^{\mu}_{\rho}\ T^{\alpha\nu}_{\verb|
|\alpha}-\delta^{\nu}_{\rho}\ T^{\alpha\mu}_{\verb| |\alpha}\right)\,.$ (5)

The modified teleparallel action for $f(T)$ gravity is given by

$I=\int d^{4}x|e|f(T)\,,$ (6)

where $|e|=\det\left(e^{A}_{\mu}\right)=\sqrt{-g}$ and the units has been
chosen so that $c=16\pi G=1$. Note that in action (6), we have overlooked any
matter contribution in the action. Varying the action in Eq. (6) with respect
to the vierbein vector field $e_{A}^{\mu}$, we obtain the equation [15]

$\frac{1}{e}\partial_{\mu}\left(eS_{A}^{\verb|
|\mu\nu}\right)f_{T}-e_{A}^{\lambda}T^{\rho}_{\verb|
|\mu\lambda}S_{\rho}^{\verb| |\nu\mu}f_{T}+S_{A}^{\verb|
|\mu\nu}\partial_{\mu}\left(T\right)f_{{}_{TT}}+\frac{1}{4}e_{A}^{\nu}f=0\,,$
(7)

where a subscript $T$ denotes differentiation with respect to $T$. We assume
the four-dimensional flat Friedmann-Lemaître-Robertson-Walker (FLRW) space-
time with the metric,

$ds^{2}=h_{\alpha\beta}dx^{\alpha}dx^{\beta}+\tilde{r}^{2}d\Omega^{2}\,,$ (8)

where $\tilde{r}=a(t)r$, $x^{0}=t$ and $x^{1}=r$ with the two-dimensional
metric $h_{\alpha\beta}={\rm diag}(-1,a^{2}(t))$. Here, $a(t)$ is the scale
factor and $d\Omega^{2}$ is the metric of two-dimensional sphere with unit
radius. In this space-time, $g_{\mu\nu}=\mathrm{diag}(-1,a^{2},a^{2},a^{2})$
and the tetrad components $e^{A}_{\mu}=(1,a,a,a)$ yield the exact value of
torsion scalar

$T=-6H^{2},$ (9)

where $H=\dot{a}/a$ is the Hubble parameter and the dot denotes the time
derivative of $\partial/\partial t$.

In the flat FLRW background, it follows from Eq. (7) that the modified
Friedmann equations are given by [15]

$\displaystyle 12f_{{}_{T}}H^{2}+f=0\,,$ (10)
$\displaystyle\dot{H}=\frac{1}{4T\,f_{{}_{TT}}+2f_{{}_{T}}}\left(-T\,f_{{}_{T}}+\frac{f}{2}\right)\,.$
(11)

It is known that $f(T)$ gravity has first-order gravitational field equation
in derivatives. Similar to $f(R)$ gravity in general relativity where the
gravitational field equation is fourth-order in derivatives, it is important
to investigate the theoretical aspects in order to examine whether $f(T)$
gravity can be a gravitational theory like $f(R)$ gravity.

In order to derive the cosmological equations in a FLRW metric, one can define
a canonical Lagrangian ${\cal L}={\cal L}(a,\dot{a},T,\dot{T})$, where ${\cal
Q}=\\{a,T\\}$ is the configuration space and ${\cal
TQ}=\\{a,\dot{a},T,\dot{T}\\}$ is the related tangent bundle on which ${\cal
L}$ is defined. The variable $a(t)$ and $T(t)$ are the scale factor and the
torsion scalar in the FLRW metric, respectively. One can use the method of the
Lagrange multipliers to set $T$ as a constraint of the dynamics. Selecting the
suitable Lagrange multiplier and integrating by parts, the Lagrangian ${\cal
L}$ becomes canonical. In our case, we have

$I=2\pi^{2}\int
dt\,a^{3}\left\\{f(T)-\lambda\left[T+6\left(\frac{\dot{a}^{2}}{a^{2}}\right)\right]\right\\},$
(12)

where $a$ is the scale factor scaled with respect to today’s value (so that
$a=\tilde{a}/\tilde{a}_{0}$ and $a(t_{0})=1$). This choice for $a$, makes it
dimensionless, whereas $[f]=M^{4}$. It is straightforward to show that, for
$f(T)=-T$, one obtains the usual Friedmann equations.

The variation with respect to $T$ of the action gives $\lambda=f_{T}$.
Therefore the previous action can be rewritten as

$I=2\pi^{2}\int
dt\,a^{3}\left\\{f-f_{T}\left[T+6\left(\frac{\dot{a}^{2}}{a^{2}}\right)\right]\right\\}\,,$
(13)

and then, integrating by parts, the point-like FLRW Lagrangian is obtained

${\cal L}=a^{3}\,(f-f_{T}\,T)-6\,f_{T}\,a\,\dot{a}^{2},$ (14)

which is a canonical function of two coupled fields, $T$ and $a$, both
depending on time $t$. The momenta conjugate to variables $a$ and $T$ are

$p_{a}=\frac{\partial{\cal L}}{\partial\dot{a}}=-12f_{T}a\dot{a}.$ (15)
$p_{T}=\frac{\partial{\cal L}}{\partial\dot{T}}=0.$ (16)

The total energy $E_{\cal L}$ corresponding to the $0,0$-Einstein equation is

$E_{{\cal L}}=-6\,f_{T}\,a\,\dot{a}^{2}-a^{3}\,(f-f_{T}\,T)=0.$ (17)

The equations of motion for $T$ and $a$ are

$a^{3}f_{{}_{TT}}\left(T+6\frac{\dot{a}^{2}}{a^{2}}\right)=0,$ (18)
$-6f_{T}\,H^{2}-12f_{T}\frac{\ddot{a}}{a}=3(f-f_{T}\,T)+12f_{{}_{TT}}\,\dot{T}\,H,$
(19)

respectively. Note that from the equations (17), (18) and (19) we can recover
the equation (10), (9) and (11), respectively. Here we have taken $f_{TT}\neq
0$.

Considering $T$ and $a$ as the variables, we see that $T$ coincides with the
definition of the torsion scalar in the FLRW metric (excluding the case
$f_{{}_{TT}}=0$). Geometrically, this is the Euler constraint of the dynamics.
Furthermore, as we will show below, constraints on the form of the function
$f(T)$ can be achieved by asking for the existence of Noether symmetries. On
the other hand, the existence of the Noether symmetries guarantees the
reduction of dynamics and the eventual solvability of the system.

##  3 The Noether symmetry

Solutions for the dynamics given by the point-like Lagrangian (14) can be
obtained by selecting cyclic variables related to some Noether symmetry [26,
27]. In principle, as a physical criterion, this approach allows one to select
$f(T)$ gravity models which are compatible with the Noether symmetry. Let
$\mathcal{L}(q^{i},\dot{q}^{i})$ be a canonical, non degenerate point-like
Lagrangian so that

$\frac{\partial\mathcal{L}}{\partial U}=0,\hskip
42.67912ptdetH_{ij}\equiv\left\|\frac{\partial^{2}\mathcal{L}}{\partial\dot{q}^{i}\partial\dot{q}^{j}}\right\|\neq
0,$ (20)

where $H_{ij}$ is the Hessian matrix related to the Lagrangian $\mathcal{L}$
and a dot denotes derivative with respect to the affine parameter $U$, namely
the cosmic time t. The Lagrangian in analytic mechanics is of the form

$\mathcal{L}=T(\textbf{q},\dot{\textbf{q}})-V(\textbf{q}),$ (21)

where T and V are the positive definite quadratic ‘kinetic energy’ and
‘potential energy’, respectively. The energy function associated with
$\mathcal{L}$ is:

$E_{\mathcal{L}}\equiv\frac{\partial\mathcal{L}}{\partial\dot{q}^{i}}-\mathcal{L},$
(22)

which is the total energy $T+V$ and the constant of motion. Since our
cosmological problem has a finite number of degrees of freedom, we are going
to consider only point transformations. Any invertible transformation of the
‘generalized positions’ $Q^{i}=Q^{i}(\textbf{q})$ induces a a local
transformation of the ‘generalized velocities’ such that:

$\dot{Q}^{i}(\textbf{q})=\frac{\partial Q^{i}}{\partial q^{j}}\dot{q}^{j};$
(23)

the matrix $\mathcal{J}=\left\|\partial Q^{i}/\partial q^{j}\right\|$ is the
Jacobian of the transformation on the positions, and it is assumed to be non-
zero. A point transformation $Q^{i}=Q^{i}(\textbf{q})$ can depend on a (or
more than one) parameter. As starting point, we can assume that a point
transformation depends on a parameter $\epsilon$, so that
$Q^{i}=Q^{i}(\textbf{q},\epsilon)$, and that it gives rise to a one-parameter
Lie group. For infinitesimal values of $\epsilon$, the transformation is then
generated by a vector field: for instance, $\partial/\partial x$ is a
translation along the $x$ axis while $x(\partial/\partial
y)-y(\partial/\partial x)$ is a rotation around the $z$ axis and so on.

In general, an infinitesimal point transformation is represented by a generic
vector field on $Q$

$\textbf{X}=\alpha^{i}(\textbf{q})\frac{\partial}{\partial q^{i}}.$ (24)

The induced transformation (23) is then represented by

$\textbf{X}^{c}=\alpha^{i}\frac{\partial}{\partial
q^{i}}+\left(\frac{d}{d\lambda}\alpha^{i}(\textbf{q})\right)\frac{\partial}{\partial\dot{q}^{j}}.$
(25)

A function $F(\textbf{q},\dot{\textbf{q}})$ is invariant under the
transformation X if

$L_{X}F\equiv\alpha^{i}(\textbf{q})\frac{\partial F}{\partial
q^{i}}+\left(\frac{d}{d\lambda}\alpha^{i}(\textbf{q})\right)\frac{\partial}{\partial\dot{q}^{j}}F=0,$
(26)

where $L_{X}F$ is the Lie derivative of F. If $L_{X}\mathcal{L}=0$, X is then
a symmetry for the dynamics derived by $\mathcal{L}$. Let us now consider a
Lagrangian $\mathcal{L}$ and its Euler-Lagrange equations

$\frac{d}{d\lambda}\frac{\partial\mathcal{L}}{\partial\dot{q}^{j}}-\frac{\partial\mathcal{L}}{\partial
q^{j}}=0.$ (27)

Contracting (27) with $\alpha^{i}$s leads to

$\alpha^{j}\left(\frac{d}{d\lambda}\frac{\partial\mathcal{L}}{\partial\dot{q}^{j}}-\frac{\partial\mathcal{L}}{\partial
q^{j}}\right)=0.$ (28)

Using the total derivative relation as

$\alpha^{j}\frac{d}{d\lambda}\frac{\partial\mathcal{L}}{\partial\dot{q}^{j}}=\frac{d}{d\lambda}\left(\alpha^{j}\frac{\partial\mathcal{L}}{\partial\dot{q}^{j}}\right)-\left(\frac{d\alpha^{j}}{d\lambda}\right)\frac{\partial\mathcal{L}}{\partial\dot{q}^{j}},$
(29)

we obtain from equation (28) that

$\frac{d}{d\lambda}\left(\alpha^{j}\frac{\partial\mathcal{L}}{\partial\dot{q}^{j}}\right)=L_{X}\mathcal{L}.$
(30)

The immediate consequence is the Noether theorem which states: if
$L_{X}\mathcal{L}=0$, then the function

$\Sigma_{0}=\alpha^{k}\frac{\partial\mathcal{L}}{\partial\dot{q}^{k}},$ (31)

is a constant of motion. Equation (31) can be expressed independently of
coordinates as a contraction of X by a Cartan 1-form

$\theta_{\mathcal{L}}\equiv\frac{\partial\mathcal{L}}{\partial\dot{q}^{j}}dq^{j}.$
(32)

For a generic vector field $\textbf{Y}=y^{i}\partial/\partial x^{i}$, and
1-form $\beta=\beta_{i}dx^{i}$, we have by definition
$i_{Y}\beta=y^{i}\beta_{i}$. Thus equation (31) can be expressed as:

$i_{X}\theta_{\mathcal{L}}=\Sigma_{0}.$ (33)

Through a point transformation, the vector field X becomes:

$\tilde{\textbf{X}}=\left(i_{X}dQ^{k}\right)\frac{\partial}{\partial
Q^{k}}+\left(\frac{d}{d\lambda}\left(i_{X}dQ^{k}\right)\right)\frac{\partial}{\partial\dot{Q}_{k}}.$
(34)

If X is a symmetry and we choose a point transformations such that

$i_{X}dQ^{1}=1;\hskip 34.14322pti_{X}dQ^{i}=0\hskip 28.45274pti\neq 1,$ (35)

we get

$\tilde{\textbf{X}}=\frac{\partial}{\partial Q^{1}};\hskip
34.14322pt\frac{\partial\mathcal{L}}{\partial Q^{1}}=0.$ (36)

Consequently, $Q^{1}$ is a cyclic coordinate, related to conserved quantities,
reducing the dynamics of the system to a manageable one. Furthermore, the
change of coordinates given by (35) is not unique. In general, the solution of
equation (35) is not defined on the whole space, rather, it is local in the
sense explained above.

Following [24, 25], we define the Noether symmetry induced on the model by a
vector field $X$ on the tangent space $T{\cal
Q}=\left(a,T,\dot{a},\dot{T}\right)$ of the configuration space ${\cal
Q}=\left(a,T\right)$ of Lagrangian (14)

$X=\alpha\frac{\partial}{\partial a}+\beta\frac{\partial}{\partial
T}+\frac{d\alpha}{dt}\frac{\partial}{\partial\dot{a}}+\frac{d\beta}{dt}\frac{\partial}{\partial\dot{T}},$
(37)

such that the Lie derivative of the Lagrangian with respect to this vector
field vanishes

$L_{X}{\cal L}=0.$ (38)

In (37), $\alpha$ and $\beta$ are functions of $a$ and $T$ and $\frac{d}{dt}$
represents the Lie derivative along the dynamical vector field, that is,

$\frac{d}{dt}=\dot{a}\frac{\partial}{\partial
a}+\dot{T}\frac{\partial}{\partial T}.$ (39)

It is easy to find the constants of motion corresponding to such a symmetry.
Indeed, equation (38) can be rewritten as

$L_{X}{\cal L}=\left(\alpha\frac{\partial{\cal L}}{\partial
a}+\frac{d\alpha}{dt}\frac{\partial{\cal
L}}{\partial\dot{a}}\right)+\left(\beta\frac{\partial{\cal L}}{\partial
T}+\frac{d\beta}{dt}\frac{\partial{\cal L}}{\partial\dot{T}}\right)=0.$ (40)

Noting that $\frac{\partial{\cal L}}{\partial q}=\frac{dp_{q}}{dt}$, we have

$\left(\alpha\frac{dp_{a}}{dt}+\frac{d\alpha}{dt}p_{a}\right)+\left(\beta\frac{dp_{T}}{dt}+\frac{d\beta}{dt}p_{T}\right)=0,$
(41)

which yields

$\frac{d}{dt}\left(\alpha p_{a}+\beta p_{T}\right)=0.$ (42)

Thus, the constants of motion are

$Q=\alpha p_{a}+\beta p_{T},$ (43)

whereas in $f(T)$ theory we have $p_{T}=0$ and so the corresponding constant
of motion becomes $Q=\alpha p_{a}$.

In order to obtain the functions $\alpha$ and $\beta$ we use equation (40). In
general, this equation gives a polynomial in terms of $\dot{a}^{2}$ and
$\dot{a}\dot{T}$ with coefficients being partial derivatives of $\alpha$ and
$\beta$ with respect to the configuration variables $a$ and $T$. Thus, the
resulting expression is identically equal to zero if and only if these
coefficients are zero. This leads to a system of partial differential
equations for $\alpha$ and $\beta$.

##  4 Noether symmetries in $f(T)$ cosmology

For the existence of a symmetry, we can write the following system of
equations (linear in $\alpha$ and $\beta$),

$f_{T}(\alpha+2a\,\partial_{{}_{a}}\alpha)+a\,f_{TT}\beta=0,$ (44)
$a\,f_{T}\,\partial_{{}_{T}}\alpha=0.$ (45)

which are obtained by setting the coefficients of the terms $\dot{a}^{2}$ and
$\dot{a}\dot{T}$ in $L_{\bf X}{\cal L}=0$ to zero. In order to make $L_{\bf
X}{\cal L}=0$ vanish we will also look for those particular $f$’s which, given
the Euler dynamics, also satisfy the constraint

$3\alpha\,(f-T\,f_{T})-a\,\beta\,T\,f_{TT}=0\,.$ (46)

This procedure is different from the usual Noether symmetry approach, in the
sense that now $L_{\bf X}{\cal L}=0$ will be solved not for all dynamics
(which solve the Euler-Lagrange equations), but only for those $f$ which
allows Euler solutions to solve also the constraint (46). Imposing such a
constraint on the form of $f$ will turn out to be, as we will show, a
sufficient condition to find solutions of the Euler-Lagrange equation which
also possess a constant of motion, i.e. a Noether charge. As we shall see
later, the system (44) and (45) can be solved exactly. Having a non-trivial
solution for $\alpha$ and $\beta$ for this system, one finds a constant of
motion if also the constraint (46) is satisfied. A solution of (44) and (45)
exists if explicit forms of $\alpha$, $\beta$ are found. If, at least one of
them is different from zero, a Noether symmetry exists.

When $f_{T}\neq 0$, from equation (45) we get

$\partial_{{}_{T}}\alpha=0\,\Rightarrow\alpha=\alpha(a).$ (47)

Thus, from this equation it can be seen $\alpha$ only depend on $a$. On the
other hand, we can obtain $\alpha$ from equation (46) as follows

$\alpha(a)=\frac{a\,f_{{}_{TT}}T}{3(f-T\,f_{{}_{T}})}\beta(a,T),$ (48)

If one uses this expression in equation (44) to eliminate $\alpha(a)$, one
obtains

$\frac{\partial\beta}{\partial a}=\frac{-3f}{2a\,T\,f_{{}_{T}}}\beta(a,T).$
(49)

To solve this equation we assume that the function $\beta(a,T)$ can be written
in the form $\beta(a,T)=A(a)B(T)$, where $A$ and $B$ are separate functions of
$a$ and $T$, respectively. Substituting this ansatz for $\beta(a,T)$ into
equation (49), we obtain

$\frac{2a}{A}\frac{dA}{da}=-\frac{3f}{T\,f_{{}_{T}}}.$ (50)

Since the left-hand side of this equation is a function of $a$ only while the
right-hand side is a function of $T$, we should have

$-\frac{3f}{T\,f_{{}_{T}}}=C=\mbox{Const},$ (51)

which results in

$f(T)=f_{{}_{0}}T^{(-\frac{3}{C})}.$ (52)

On the other hand equation (50), with its right hand-side equal to $C$, has
the solution

$\frac{2a}{A}\frac{dA}{da}=C\Rightarrow A(a)=a^{C/2}.$ (53)

Now, using (52) and $\beta(a,T)=a^{C/2}B(T)$ in equation (48), we find

$\alpha(a)=\frac{1}{C}T^{-1}\,a^{\frac{C}{2}+1}B(T).$ (54)

Since $\alpha(a)$ should be a function of $a$ only, from the above expression
for $\alpha(a)$ we can write

$B(T)=T,$ (55)

and thus we obtain

$\beta(a,T)=T\,a^{C/2},\hskip
14.22636pt\alpha(a)=\frac{1}{C}a^{\frac{C}{2}+1}.$ (56)

To obtain the corresponding cosmology resulting from this type of $f(T)$, we
note that the existence of Noether symmetry implies the existence of a
constant of motion $Q=\alpha p_{a}$. Hence, using equation (15) we have

$Q=36\frac{a^{(\frac{C}{2}+2)}\dot{a}\,T^{-(1+\frac{3}{C})}}{C}.$ (57)

On the other hand, by taking $f_{{}_{TT}}\neq 0$ from equation (18) we obtain

$T=-6\frac{\dot{a}^{2}}{a^{2}}.$ (58)

By using equations (57) and (58) we can obtain scale factor as following

$a(t)\sim t^{\frac{-2}{C}}.$ (59)

Therefore, in the context of $f(T)=f_{{}_{0}}T^{-(\frac{3}{C})}$ cosmology,
the universe evolves with a power law expansion. It is seen from (59) that the
condition under which the universe expands is $C<0$. The deceleration
parameter as a function of $C$ is therefore given by

$q(C)=-\left(1+\frac{C}{2}\right)$ (60)

The condition for acceleration is $q(C)<0$, thus we have $C>-2$. As it can be
seen, for $C\rightarrow 0$ we have $q\rightarrow-1$, that is the universe
finally approaches the eternal de Sitter phase with infinite acceleration,
however, the accelerated expansion occurs in $-2<C<0$.

##  5 Dark energy equation of state and age of the Universe in $f(T)$

In this section, we consider the accelerated expansion of the universe in the
context of $f(T)$ theory, without introducing the mysteries fluid the so-
called “dark energy” with a negative equation of state (EOS) parameter $w$.
Let us then start from equation (59) and find the Hubble parameter as a
function of the redshift $z$ as

$H(z)=-\frac{2}{C}H_{0}(1+z)^{-\frac{C}{2}},$ (61)

where $a_{0}/a=1+z$ with $a_{0}$ and $H_{0}$ being the values of the parameter
at the present epoch. It is worth noting that equation (61) is the same as
that derived from the standard Friedmann equation with $w=-2/3$. Now, we may
write the Friedmann equation in a formal fashion which would encapsulate any
modification to the standard Friedmann equation in the last term regardless of
its nature [28] that is

$H^{2}/H^{2}_{0}=\Omega_{m}(1+z)^{3}+\delta H^{2}/H^{2}_{0},$ (62)

where $\Omega_{m}=\rho/\rho_{0c}$, $\rho_{0c}=3H^{2}_{0}$. Also, defining the
effective EOS, denoted by $w_{{}_{\rm eff}}(z)$, as

$w_{{}_{\rm eff}}(z)=-1+\frac{1}{3}\frac{d\ln\delta H^{2}}{d\ln(1+z)},$ (63)

we can calculate $w_{{}_{\rm eff}}(z)$ using equations (61), (62) and (63)
with the result

$w_{{}_{\rm
eff}}=-1+\frac{1}{3}\frac{-\frac{1}{C}(1+z)^{-c}-3\Omega_{m}(1+z)^{3}}{\frac{4}{C^{2}}(1+z)^{-c}-\Omega_{m}(1+z)^{3}}.$
(64)

Figure 1 shows the behavior of the effective EOS parameter, $w_{{}_{\rm
eff}}$, as a function of $z$. As it can be seen, for $C=-1$ ($-2<C\leq-1$) and
$\Omega_{m}=0.33$ we have $w_{{}_{\rm eff}}\lessapprox-1$, which is the
characteristic of one type of dark energy, the so-called phantom and from
equation (64), for $C\rightarrow 0^{-}$, we have $w_{{}_{\rm
eff}}\rightarrow-1$.

Figure 1: The behavior of effective EOS, $w_{{}_{\rm eff}}$, as a function of
$z$. When $C=-1$ and $\Omega_{m}=0.33$, we get a phantom behavior in the
model.

To continue we consider the age of the Universe in $f(T)$ model. Thus, the age
of the matter dominated Universe in FLRW models is given by

$H_{0}t_{0}=\frac{2}{3}\frac{1}{\sqrt{1-\Omega_{m}}}\ln\left[\frac{1+\sqrt{1-\Omega_{m}}}{\sqrt{\Omega_{m}}}\right]$
(65)

where $H_{0}^{-1}=9.8\times 10^{9}h^{-1}$ years and the dimensionless
parameter $h$, according to resent data, is about $0.7$. Hence, in the flat
matter dominated universe with $\Omega_{total}=1$ the age of the universe
would be only 9.3 Gyr, whereas the oldest globular clusters yield an age of
about $\sim 13.5$ Gyr [29]. We obtain the age of the universe for our model
from equations (10), (52) and (59) by taking matter in the Friedmann equations
as follows

$H_{0}t_{0}=\frac{-2}{C}\frac{1}{\sqrt{1-\Omega_{m}}}\ln\left[\frac{1+\sqrt{1-\Omega_{m}}}{\sqrt{\Omega_{m}}}\right]$
(66)

For a flat, matter dominated Universe with $\Omega_{m}\approx 0.3$ and $C=-3$
we have a prediction for the age of the Universe of about $13.2$ Gyr. It seems
that the age of the universe in our model is longer than the FLRW model.
Figure 2 shows the behavior of the dimensionless age parameter, $H_{0}t_{0}$,
as a function of $\Omega_{m}$ for different values of $C$.

Figure 2: $H_{0}t_{0}$ as a function of $\Omega_{m}$ for $C=-1.9$ (solid
line), $C=-3$ (dashed line) and $C=-6$ (dotted line). Figure shows that for a
fixed value of $\Omega_{m}$ the predicted age of the universe is longer for
larger values of $C$.

Note added

Before sending our present work to the Archive (on 13 December 2011), we
became aware of the fact that there already existed a paper arXiv:1112.2270
[gr-qc][30] (sent to the Archive on 10 December 2011), in which the authors
had addressed similar issues. Although the content and results of the present
work overlap with the ones in [30], it is worth noting that our work was
completed independently.

## Acknowledgment

This work has been supported financially by Research Institute for Astronomy
and Astrophysics of Maragha (RIAAM) under research project No. 1/2361.

##  6 Conclusions

In this paper we have studied a generic $f(T)$ cosmological model by using
Noether symmetry approach. We have taken the background geometry as a flat
FLRW metric and derived the general equations of motion in this background.
The phase space has been constructed by taking the scale factor $a$ and
torsion scalar $T$ as two independent dynamical variables. The Lagrangian of
the model in the configuration space spanned by $\left\\{a,T\right\\}$ is so
constructed that its variation with respect to these dynamical variables
yields the correct field equations. The existence of Noether symmetry implies
that the Lie derivative of this Lagrangian with respect to the infinitesimal
generator of the desired symmetry vanishes. By applying this condition to the
Lagrangian of the model, we have obtained the explicit form of the
corresponding $f(T)$ function. We have shown that this form of $f(T)$ results
in a power law expansion for the scale factor and the accelerated expansion
occurs in $-2<C<0$

## References

  * [1] S. Perlmutter et al. [SNCP Collaboration], Astrophys. J. 517 (1999) 565.
  * [2] A. G. Riess et al. [SNST Collaboration], Astron. J. 116 (1998) 1009 .
  * [3] D. N. Spergel et al. [WMAP Collaboration], Astrophys. J. Suppl. 148 (2003) 175 .
  * [4] E. Komatsu et al. [WMAP Collaboration], Astrophys. J. Suppl. 192 (2011) 18 .
  * [5] M. Tegmark et al., Phys. Rev. D 69 (2004) 103501.
  * [6] D. J. Eisenstein et al., Astrophys. J. 633 (2005) 560 .
  * [7] B. Jain and A. Taylor, Phys. Rev. Lett. 91 (2003) 141302 .
  * [8] E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15 (2006) 1753 .
  * [9] M. Li, X. D. Li, S. Wang and Y. Wang, Commun. Theor. Phys. 56 (2011) 525, arXiv:1103.5870.
  * [10] A. De Felice and S. Tsujikawa, Living Rev. Rel. 13 (2010) 3 .
  * [11] T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, arXiv:1106.2476.
  * [12] R. Weitzenböck, Invarianten Theorie, (Nordhoff, Groningen, 1923).
  * [13] F. W. Hehl, P. Von Der Heyde, G. D. Kerlick and J. M. Nester, Rev. Mod. Phys. 48 (1976) 393.
  * [14] A. Einstein, Sitzungsber. Preuss. Akad. Wiss. Phys. Math. Kl., 217 (1928); 401 (1930).
  * [15] G. R. Bengochea and R. Ferraro, Phys. Rev. D 79 (2009) 124019.
  * [16] R. Ferraro and F. Fiorini, Phys. Rev. D 75 (2007) 084031.
  * [17] S. H. Chen, J. B. Dent, S. Dutta and E. N. Saridakis, Phys. Rev. D 83, (2011) 023508.
  * [18] P. Wu and H. W. Yu, Eur. Phys. J. C 71 (2011) 1552.
  * [19] J. B. Dent, S. Dutta and E. N. Saridakis, JCAP 1101 (2011) 009.
  * [20] T. Wang, Phys. Rev. D 84 (2011) 024042.
  * [21] B. Li, T. P. Sotiriou and J. D. Barrow, Phys. Rev. D 83 (2011) 064035.
  * [22] R. Aldrovandi and J. G. Pereira, An Introduction to Teleparallel Gravity. unpublished coursenotes. www.ift.unesp.br/gcg/tele.pdf
  * [23] J. C. C. de Souza and V. Faraoni, Class. Quantum Grav. 24 (2007) 3637, arXiv: 0706.1223.
  * [24] M. Demianski, R. de Ritis, C. Rubano and P. Scudellaro, Phys. Rev. D 46 (1992) 1391\.
  * [25] B. Vakili, Phys. Lett. B 669 (2008) 206, arXiv:0809.4591.
  * [26] S. Capozziello, A. De Felice, JCAP 0808 (2008) 016.
  * [27] B. Vakili, N. Khosravi and H. R. Sepangi, Class. Quantum Grav. 24 (2007) 931, arXiv: gr-qc/0701075.
  * [28] E. V. Linder and A. Jenkins, Mon. Not. R. Astron. Soc. 346 (2003) 573, arXiv: astroph/ 0305286\.
  * [29] L. M. Krauss and B. Chaboyer, Science 299 (2003) 65.
  * [30] H. Wei, X-J. Guo, L-F. Wang, Phys. Lett. B 707 (2012) 298, arXiv:1112.2270.