Kuan Liu

Aurélien Bellet

Fei Sha

Similarity Learning for High-Dimensional Sparse Data

A good measure of similarity between data points is crucial to many tasks in machine learning. Similarity and metric learning methods learn such measures automatically from data, but they do not scale well respect to the dimensionality of the data. In this paper, we propose a method that can learn efficiently similarity measure from highdimensional sparse data. The core idea is to parameterize the similarity measure as a convex combination of rank-one matrices with specific sparsity structures. The parameters are then optimized with an approximate Frank-Wolfe procedure to maximally satisfy relative similarity constraints on the training data. Our algorithm greedily incorporates one pair of features at a time into the similarity measure, providing an efficient way to control the number of active features and thus reduce overfitting. It enjoys very appealing convergence guarantees and its time and memory complexity depends on the sparsity of the data instead of the dimension of the feature space. Our experiments on realworld high-dimensional datasets demonstrate its potential for classification, dimensionality reduction and data exploration.

INTRODUCTION

In many applications, such as text processing, computer vision or biology, data is represented as very high-dimensional but sparse vectors. The ability to compute meaningful similarity scores between these objects is crucial to many tasks, such as classification, clustering or ranking. However, handcrafting a relevant similarity measure for such data is challenging because it is usually the case that only a small, often unknown subset of features is actually relevant to the task at hand. For instance, in drug discovery, chemical compounds can be represented as sparse features describing their 3D properties, and only a few of them play an role in determining whether the compound will bind to a target receptor [START_REF] Guyon | Result Analysis of the NIPS 2003 Feature Selection Challenge[END_REF]. In text classification, where each document is represented as a sparse bag of words, only a small subset of the words is generally sufficient to discriminate among documents of different topics.

A principled way to obtain a similarity measure tailored to the problem of interest is to learn it from data. This line of research, known as similarity and distance metric learning, has been successfully applied to many application domains (see [START_REF] Kulis | Metric Learning: A Survey[END_REF]Bellet et al., 2013, for recent surveys). The basic idea is to learn the parameters of a similarity (or distance) function such that it satisfies proximity-based constraints, requiring for instance that some data instance x be more similar to y than to z according to the learned function. However, similarity learning typically requires estimating a matrix with O(d 2 ) entries (where d is the data dimension) to account for correlation between pairs of features. For high-dimensional data (say, d > 10 4 ), this is problematic for at least three reasons: (i) training the metric is computationally expensive (quadratic or cubic in d), (ii) the matrix may not even fit in memory, and (iii) learning so many parameters is likely to lead to severe overfitting, especially for sparse data where some features are rarely observed.

To overcome this difficulty, a common practice is to first project data into a low-dimensional space (using PCA or random projections), and then learn a similarity function in the reduced space. Note that the projection intertwines useful features and irrelevant/noisy ones. Moreover, it is also difficult to interpret the reduced feature space, when we are interested in discovering what features are more important than others for discrimination.

In this paper, we propose a novel method to learn a bilinear similarity function S M (x, x ) = x T M x directly in the original high-dimensional space while avoiding the above-mentioned pitfalls. The main idea combines three ingredients: the sparsity of the data, the parameterization of M as a convex combination of rank-one matrices with special sparsity structures, and an approximate Frank-Wolfe procedure [START_REF] Frank | An algorithm for quadratic programming[END_REF][START_REF] Kenneth | Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm[END_REF][START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF] to learn the similarity parameters. The resulting algorithm iteratively and greedily incorporates one pair of features at a time into the learned similarity, providing an efficient way to ignore irrelevant features as well as to guard against overfitting through early stopping. Our method has appealing approximation error guarantees, time and memory complexity independent of d and outputs extremely sparse similarity functions that are fast to compute and to interpret.

The usefulness of the proposed approach is evaluated on several datasets with up to 100,000 features, some of which have a large proportion of irrelevant features. To the best of our knowledge, this is the first time that a full similarity or distance metric is learned directly on such high-dimensional datasets without first reducing dimensionality. Our approach significantly outperforms both a diagonal similarity learned in the original space and a full similarity learned in a reduced space (after PCA or random projections). Furthermore, our similarity functions are extremely sparse (in the order of 0.0001% of nonzero entries), using a sparse subset of features and thus providing more economical analysis of the resulting model (for example, examining the importance of the original features and their pairwise interactions).

The rest of this paper is organized as follows. Section 2 briefly reviews some related work. Our approach is described in Section 3. We present our experimental results in Section 4 and conclude in Section 5.

RELATED WORK

Learning similarity and distance metric has attracted a lot of interests. In this section, we review previous efforts that focus on efficient algorithms for highdimensional data -a comprehensive survey of existing approaches can be found in [START_REF] Bellet | A Survey on Metric Learning for Feature Vectors and Structured Data[END_REF].

A majority of learning similarity has focused on learning either a Mahalanobis distance d M (x, x ) = (xx ) T M (x -x ) where M is a symmetric positive semi-definite (PSD) matrix, or a bilinear similarity S M (x, x ) = x T M x where M is an arbitrary d × d matrix. In both cases, it requires estimating O(d 2 ) parameters, which is undesirable in the high-dimensional setting. Virtually all existing methods thus resort to dimensionality reduction (such as PCA or random projections) to preprocess the data when it has more than a few hundred dimensions, thereby incurring a po-tential loss of performance and interpretability of the resulting function (see e.g., [START_REF] Davis | Information-theoretic metric learning[END_REF][START_REF] Kilian | Distance Metric Learning for Large Margin Nearest Neighbor Classification[END_REF][START_REF] Guillaumin | Is that you? Metric learning approaches for face identification[END_REF][START_REF] Ying | Distance Metric Learning with Eigenvalue Optimization[END_REF][START_REF] Wang | Parametric Local Metric Learning for Nearest Neighbor Classification[END_REF][START_REF] Lim | Robust Structural Metric Learning[END_REF][START_REF] Qi Qian | An Integrated Framework for High Dimensional Distance Metric Learning and Its Application to Fine-Grained Visual Categorization[END_REF].

There have been a few solutions to this essential limitation. The most drastic strategy is to learn a diagonal matrix M [START_REF] Schultz | Learning a Distance Metric from Relative Comparisons[END_REF][START_REF] Gao | SOML: Sparse Online Metric Learning with Application to Image Retrieval[END_REF], which is very restrictive as it amounts to a simple weighting of the features. Instead, some approaches assume an explicit low-rank decomposition M = L T L and learn L ∈ R r×d in order to reduce the number of parameters to learn [START_REF] Goldberger | Neighbourhood Components Analysis[END_REF][START_REF] Kilian | Distance Metric Learning for Large Margin Nearest Neighbor Classification[END_REF][START_REF] Kedem | Non-linear Metric Learning[END_REF]. But this results in nonconvex formulations with many bad local optima [START_REF] Kulis | Metric Learning: A Survey[END_REF] and requires to tune r carefully. Moreover, the training complexity still depends on d and can thus remain quite large. Another direction is to learn M as a combination of rank-one matrices. In [START_REF] Shen | Positive Semidefinite Metric Learning Using Boosting-like Algorithms[END_REF], the combining elements are selected greedily in a boosting manner but each iteration has an O(d 2 ) complexity. To go around this limitation, [START_REF] Shi | Sparse Compositional Metric Learning[END_REF] generate a set of rankone matrices before training and learn a sparse combination. However, as the dimension increases, a larger dictionary is needed and can be expensive to generate. Some work have also gone into sparse and/or low-rank regularization to reduce overfitting in high dimensions [START_REF] Rosales | Learning Sparse Metrics via Linear Programming[END_REF][START_REF] Qi | An Efficient Sparse Metric Learning in High-Dimensional Space via l1-Penalized Log-Determinant Regularization[END_REF][START_REF] Ying | Sparse Metric Learning via Smooth Optimization[END_REF] but those do not reduce the training complexity of the algorithm.

To the best of our knowledge, DML-eig [START_REF] Ying | Distance Metric Learning with Eigenvalue Optimization[END_REF] and its extension DML-ρ [START_REF] Cao | Distance Metric Learning Revisited[END_REF] are the only prior attempts to use a Frank-Wolfe procedure for metric or similarity learning. However, their formulation requires finding the largest eigenvector of the gradient matrix at each iteration, which scales in O(d 2 ) and is thus unsuitable for the high-dimensional setting we consider in this work.

PROPOSED APPROACH

This section introduces hdsl (High-Dimensional Similarity Learning), the approach proposed in this paper. We first describe our problem formulation (Section 3.1), then derive an efficient algorithm to solve it (Section 3.2).

Problem Formulation

In this work, we propose to learn a similarity function for high-dimensional sparse data. We assume the data points lie in some space X ⊆ R d , where d is large (d > 10 4 ), and are D-sparse on average (D d). Namely, the number of nonzero entries is typically much smaller than d. We focus on learning a similarity function S

M : X × X → R of the form S M (x, x ) = x T M x ,
where M ∈ R d×d . Note that for any M , S M can be computed in O(D 2 ) time on average.

Feasible domain Our goal is to derive an algorithm to learn a very sparse M with time and memory requirements that depend on D but not on d. To this end, given a scale parameter λ > 0, we will parameterize M as a convex combination of 4-sparse d × d bases:

M ∈ D λ = conv(B λ ), with B λ = ij P (ij) λ , N (ij) λ ,
where for any pair of features i, j ∈ {1, . . . , d}, i = j,

P (ij) λ = λ(e i + e j )(e i + e j ) T =   • • • • • • λ • λ • • • • • • • λ • λ • • • • • •   , N (ij) λ = λ(e i -e j )(e i -e j ) T =   • • • • • • λ • -λ • • • • • • • -λ • λ • • • • • •   .
The use of such sparse matrices was suggested by [START_REF] Jaggi | Sparse Convex Optimization Methods for Machine Learning[END_REF]. Besides the fact that they are instrumental to the efficiency of our algorithm (see Section 3.2), we give some additional motivation for their use in the context of similarity learning.

First, any M ∈ D λ is a convex combination of symmetric PSD matrices and is thus also symmetric PSD. Unlike many metric learning algorithms, we thus avoid the O(d 3 ) cost of projecting onto the PSD cone. Furthermore, constraining M to be symmetric PSD provides useful regularization to prevent overfitting [START_REF] Chechik | An online algorithm for large scale image similarity learning[END_REF] and allows the use of the square root of M to project the data into a new space where the dot product is equivalent to S M . Because the bases in B λ are rank-one, the dimensionality of this projection space is at most the number of bases composing M .

Second, each basis operates on two features only. In particular, S P (ij) λ (x, x ) = λ(x i x i + x j x j + x i x j + x j x i ) assigns a higher score when feature i appears jointly in x and x (likewise for j), as well as when feature i in x and feature j in y (and vice versa) co-occur. Conversely, S N (ij) λ penalizes the co-occurrence of features i and j. This will allow us to easily control the number of active features and learn a very compact similarity representation.

Finally, notice that in the context of text data represented as bags-of-words (or other count data), the bases in B λ are quite natural: they can be intuitively thought of as encoding the fact that a term i or j present in both documents makes them more similar, and that two terms i and j are associated with the same/different class or topic.

Optimization problem We now describe the optimization problem to learn the similarity parameters. Following previous work (see for instance [START_REF] Schultz | Learning a Distance Metric from Relative Comparisons[END_REF][START_REF] Kilian | Distance Metric Learning for Large Margin Nearest Neighbor Classification[END_REF][START_REF] Chechik | An online algorithm for large scale image similarity learning[END_REF], our training data consist of side information in the form of triplet constraints: T = {x t should be more similar to y t than to z t } T t=1 .

Such constraints can be built from a labeled training sample, provided directly by a domain expert, or obtained through implicit feedback such as clicks on search engine results. For notational convenience, we write A t = x t (y t -z t ) T ∈ R d×d for each constraint t = 1, . . . , T . We want to define an objective function that applies a penalty when a constraint t is not satisfied with margin at least 1, i.e. whenever A t , M = S M (x t , y t ) -S M (x t , z t ) < 1. To this end, we use the smoothed hinge loss : R → R + :

A t , M =    0 if A t , M ≥ 1 1 2 -A t , M if A t , M ≤ 0 1 2 (1 -A t , M ) 2 otherwise
, where •, • denotes the Frobenius inner product.1 Given λ > 0, our similarity learning formulation aims at finding the matrix M ∈ D λ that minimizes the average penalty over the triplet constraints in T : min

M ∈D λ f (M ) = 1 T T t=1 A t , M (1)
Due to the convexity of the smoothed hinge loss, Problem (1) involves minimizing a convex function over the convex domain D λ . In the next section, we propose a greedy algorithm to solve this problem.

Algorithm

Exact Frank-Wolfe Algorithm

We propose to use a Frank-Wolfe (FW) algorithm [START_REF] Frank | An algorithm for quadratic programming[END_REF][START_REF] Kenneth | Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm[END_REF][START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF] to learn the similarity. FW is a general procedure to minimize a convex and continuously differentiable function over a compact and convex set. At each iteration, it moves towards a feasible point that minimizes a linearization of the objective function at the current Algorithm 1 Frank Wolfe algorithm for problem (1)

1: initialize M (0) to an arbitrary B ∈ B λ 2: for k = 0, 1, 2, . . . do 3: let B (k) F ∈ arg min B∈B λ B, ∇f (M (k) ) and D (k) F = B (k) F -M (k) // compute forward direction 4: let B (k) A ∈ arg max B∈S (k) B, ∇f (M (k) ) and D (k) A = M (k) -B (k) A // compute away direction 5: if D (k) F , ∇f (M (k) ) ≤ D (k) A , ∇f (M (k) ) then 6: D (k) = D (k)
F and γ max = 1 // choose forward step 7: else 8:

D (k) = D (k) A and γ max = α (k) B (k) A /(1 -α (k) B (k) A ) // choose away step 9: end if 10: let γ (k) ∈ arg min γ∈[0,γmax] f (M (k) + γD (k) ) // perform line search 11: M (k+1) = M (k) + γ (k) D (k)
// update iterate towards direction 12: end for iterate. Note that a minimizer of this linear function must be at a vertex of the feasible domain. We will exploit the fact that in our formulation (1), the vertices of the feasible domain D λ are the elements of B λ and have special structure.

The FW algorithm applied to (1) and enhanced with so-called away steps [START_REF] Guélat | Some comments on Wolfe's away step[END_REF]) is described in details in Algorithm 1. During the course of the algorithm, we explicitly maintain a representation of each iterate M (k) as a convex combination of basis elements:

M (k) = B∈B λ α (k) B B, where B∈B λ α (k) B = 1, α (k) B ≥ 0.
We denote the set of active basis elements in M (k) as

S (k) = {B ∈ B λ : α (k)
B > 0}. The algorithm goes as follows. We initialize M (0) to a random basis element. Then, at each iteration, we greedily choose between moving towards a (possibly) new basis (forward step) or reducing the weight of an active one (away step). The extent of the step is determined by line search. As a result, Algorithm 1 adds only one basis (at most 2 new features) at each iteration, which provides a convenient way to control the number of active features and maintain a compact representation of M (k) in O(k) memory cost. Furthermore, away steps provide a way to reduce the importance of a potentially "bad" basis element added at an earlier iteration (or even remove it completely when γ (k) = γ max ). Note that throughout the execution of the algorithm, all iterates M (k) remain convex combinations of basis elements and are thus feasible. The following lemma shows that the iterates of Algorithm 1 converge to an optimal solution of (1) with a rate of O(1/k).

Lemma 1. Let λ > 0, M * be an optimal solution to (1) and

L = 1 T T t=1 A t 2 F . At any iteration k ≥ 1 of Algorithm 1, the iterate M (k) ∈ D λ satis- fies f (M (k) ) -f (M * ) ≤ 16Lλ 2 /(k + 2). Furthermore,
it has at most rank k + 1 with 4(k + 1) nonzero entries, and uses at most 2(k + 1) distinct features.

Proof. The result follows from the analysis of the general FW algorithm [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF], the fact that f has L-Lipschitz continuous gradient and observing that diam

• F (D λ ) = √ 8λ.
Note that the optimality gap in Lemma 1 is independent from d. This means that Algorithm 1 is able to find a good approximate solution based on a small number of features, which is very appealing in the high-dimensional setting.

Complexity Analysis

We now analyze the time and memory complexity of Algorithm 1. Observe that the gradient has the form

∇f (M ) = 1 T T t=1 G t ,
where

G t =    0 if A t , M ≥ 1 -A t if A t , M ≤ 0 ( A t , M -1) A t otherwise .
The structure of the algorithm's updates is crucial to its efficiency: since M (k+1) is a convex combination of M (k) and a 4-sparse matrix B (k) , we can efficiently compute most of the quantities of interest through careful book-keeping.

In particular, storing M (k) at iteration k requires O(k) memory. We can also recursively compute A t , M (k+1) for all constraints in only O(T ) time and O(T ) memory based on A t , M (k) and A t , B (k) . This allows us, for instance, to efficiently compute the objective value as well identify the set of satisfied constraints (those with A t , M (k) ≥ 1) and ignore them when computing the gradient. Finding the away direction at iteration k can be done in O(T k) time. For the line search, we use a bisection algorithm to find a

+ T k) Õ(T D 2 + k) Mini-batch Õ(M D 2 + T k) Õ(T + M D 2 + k) Heuristic Õ(M D + T k) Õ(T + M D + k)
root of the gradient of the 1-dimensional function of γ, which only depends on A t , M (k) and A t , B (k) , both of which are readily available. Its time complexity is O(T log 1 ) where is the precision of the linesearch, and requires constant memory.

The bottleneck is to find the forward direction. Indeed, sequentially considering each basis element is intractable as it takes O(T d2 ) time. A more efficient strategy is to sequentially consider each constraint, which requires O(T D 2 ) time and O(T D 2 ) memory.

The overall iteration complexity of Algorithm 1 is given in Table 1.

Approximate Forward Step

Finding the forward direction can be expensive when T and D are both large. We propose two strategies to alleviate this cost by finding an approximately optimal basis (see Table 1 for iteration complexity).

Mini-Batch Approximation Instead of finding the forward and away directions based on the full gradient at each iteration, we estimate it on a mini-batch of M T constraints drawn uniformly at random (without replacement). The complexity of finding the forward direction is thus reduced to O(M D 2 ) time and O(M D 2 ) memory. Under mild assumptions, concentration bounds such as Hoeffding's inequality without replacement [START_REF] Serfling | Probability inequalities for the sum in sampling without replacement[END_REF] can be used to show that with high probability, the deviation between the "utility" of any basis element B on the full set of constraints and its estimation on the mini-batch, namely:

1 M M t=1 B, G t - 1 T T t=1 B, G t , decreases as O(1/ √ M ).
In other words, the minibatch variant of Algorithm 1 finds a forward direction which is approximately optimal. The FW algorithm is known to be robust to this setting, and convergence guarantees similar to Lemma 1 can be obtained following [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF][START_REF] Freund | New Analysis and Results for the Conditional Gradient Method[END_REF].

Fast Heuristic To avoid the quadratic dependence on D, we propose to use the following heuristic to find a good forward basis. We first pick a feature i ∈ {1, . . . , d} uniformly at random, and solve the linear problem over the restricted set j {P

(ij) λ , N (ij) λ }.
We then fix j and solve the problem again over the set

EXPERIMENTS

In this section, we present experiments to study the performance of HDSL in classification, dimensionality reduction and data exploration against competing approaches.

Experimental Setup

Datasets We report experimental results on several real-world classification datasets with up to 100,000 features. Dorothea and dexter come from the NIPS 2003 feature selection challenge [START_REF] Guyon | Result Analysis of the NIPS 2003 Feature Selection Challenge[END_REF] and are respectively pharmaceutical and text data with predefined splitting into training, validation and test sets. They both contain a large proportion of noisy/irrelevant features. Reuters CV1 is a popular text classification dataset with bag-of-words representation. We use the binary version from the LIB-SVM dataset collection 2 (with 60%/20%/20% random splits) and the 4-classes version (with 40%/30%/30% random splits) introduced in [START_REF] Cai | Manifold Adaptive Experimental Design for Text Categorization[END_REF]. Detailed information on the datasets and splits is given in Table 2. All datasets are normalized such that each feature takes values in [0, 1].

Competing Methods We compare the proposed approach (hdsl) to several methods:

• identity: The standard dot product as a baseline, which corresponds to using M = I.

• diag: Diagonal similarity learning (i.e., a weighting of the features), as done in [START_REF] Gao | SOML: Sparse Online Metric Learning with Application to Image Retrieval[END_REF].

We obtain it by minimizing the same loss as in hdsl with 2 and 1 regularization, i.e., min

w∈R d f (w) = 1 T T t=1 A t , diag(w) + λΩ(w),
where Ω(w) ∈ { w 2 2 , w 1 } and λ is the regularization parameter. Optimization is done using (proximal) gradient descent. • rp+oasis: Similarity learning in random projected space. Given r d, let R ∈ R d×r be a matrix where each entry r ij is randomly drawn from N (0, 1). For each data instance x ∈ R d , we generate x = 1 √ r R T x ∈ R r and use this reduced data in OASIS [START_REF] Chechik | An online algorithm for large scale image similarity learning[END_REF], a fast online method to learn a bilinear similarity from triplet constraints.

• pca+oasis: Similarity learning in PCA space.

Same as rp+oasis, except that PCA is used instead of random projections to project the data into R r .

• svm: Support Vector Machines. We use linear SVM, which is known to perform well for sparse high-dimensional data [START_REF] Caruana | An empirical evaluation of supervised learning in high dimensions[END_REF], with 2 and 1 regularization. We also use nonlinear SVM with the polynomial kernel (2nd and 3rd degree) popular in text classification [START_REF] Chang | Training and Testing Low-degree Polynomial Data Mappings via Linear SVM[END_REF]. The SVM models are trained using liblinear [START_REF] Fan | LIBLINEAR: A Library for Large Linear Classification[END_REF] and libsvm [START_REF] Chang | LIBSVM : a library for support vector machines[END_REF] with 1vs1 paradigm for multiclass.

Training Procedure For all similarity learning algorithms, we generate 15 training constraints for each instance by identifying its 3 target neighbors (nearest neighbors with same label) and 5 impostors (nearest neighbors with different label), following [START_REF] Kilian | Distance Metric Learning for Large Margin Nearest Neighbor Classification[END_REF]. Due to the very small number of training instances in dexter, we found that better performance is achieved using 20 constraints per instance x, each of them constructed by randomly drawing a point from the class of x and a point from a different class. All parameters are tuned using the accuracy on the validation set. For hdsl, we use the fast heuristic described in Section 3.2.3 and tune the scale parameter λ ∈ {1, 10, . . . , 10 9 }. The regularization parameters of diag and svm are tuned in {10 -9 , . . . , 10 8 } and the "aggressiveness" parameter of OASIS is tuned in {10 -9 , . . . , 10 2 }.

Results

Classification Performance We first investigate the performance of each similarity learning approach in k-NN classification (k was set to 3 for all experiments). For rp+oasis and pca+oasis, we choose the dimension r of the reduced space based on the accuracy of the learned similarity on the validation set, limiting our search to r ≤ 2000 because OASIS is extremely slow beyond this point. 3 Similarly, we use the performance on validation data to do early stopping in hdsl, which also has the effect of restricting the number of features used by the learned similarity.

Table 3 shows the k-NN classification performance. First, notice that rp+oasis often performs worse than identity, which is consistent with previous observations that a large number of random projections may be needed to obtain good performance [START_REF] Fradkin | Experiments with random projections for machine learning[END_REF]. pca+oasis gives much better results, but is generally outperformed by a simple diagonal similarity learned directly in the original highdimensional space. hdsl, however, outperforms all other algorithms on these datasets, including diag. This shows the good generalization performance of the proposed approach, even though the number of training samples is sometimes very small compared to the number of features, as in dexter and dorothea. It also shows the relevance of encoding "second order" information (pairwise interactions between the original features) in the similarity instead of simply considering a feature weighting as in diag.

Table 4 shows the comparison with SVMs. Interestingly, hdsl with k-NN outperforms all SVM variants on dexter and dorothea, both of which have a large proportion of irrelevant features. This shows that its greedy strategy and early stopping mechanism achieves better feature selection and generalization than the 1 version of linear SVM. On the other two datasets, hdsl is competitive with SVM, although it is outperformed slightly by one variant (svm-poly-3 on rcv1 2 and svm-linear-2 on rcv1 4), both of which rely on all features.

Feature Selection and Sparsity We now focus on the ability of hdsl to perform feature selection and more generally to learn sparse similarity functions. To better understand the behavior of hdsl, we show in Figure 1 the number of selected features as a function of the iteration number for two of the datasets. Remember that at most two new features can be added at each iteration. Figure 1 shows that hdsl incorporates many features early on but tends to eventually converge to a modest fraction of features (the same observation holds for the other two datasets). This may explain why hdsl does not suffer much from overfitting even when training data is scarce as in dexter.

Another attractive characteristic of hdsl is its ability to learn a matrix that is sparse not only on the diagonal but also off-diagonal (the proportion of nonzero entries is in the order of 0.0001% for all datasets). In other words, it only relies on a few relevant pairwise interactions between features. Figure 2 shows two examples, where we can see that hdsl is able to exploit the product of two features as either a positive or negative contribution to the similarity score. This opens the door to an analysis of the importance of pairs of features (for instance, word co-occurrence) for the application at hand. Finally, the extreme sparsity of the matrices allows very fast similarity computation.

Finally, it is also worth noticing that hdsl uses significantly less features than diag-1 (see numbers in brackets in Table 3). We attribute this to the extra modeling capability brought by the non-diagonal similarity observed in Figure 2. 4Dimension Reduction We now investigate the potential of hdsl for dimensionality reduction. Recall that hdsl learns a sequence of PSD matrices M (k) . We can use the square root of M (k) to project the data into a new space where the dot product is equivalent to S M (k) in the original space. The dimension of the projection space is equal to the rank of M (k) , which is upper bounded by k +1. A single run of hdsl can thus be seen as incrementally building projection spaces of increasing dimensionality.

To assess the dimensionality reduction quality of hdsl (measured by k-NN classification error on the test set), we plot its performance at various iterations during the runs that gave the results in Table 3. We compare it to two standard dimensionality reduction techniques: random projection and PCA. We also evaluate rp+oasis and pca+oasis, i.e., learn a similarity with OASIS on top of the RP and PCA features.5 Note that OASIS was tuned separately for each projection size, making the comparison a bit unfair to hdsl. The results are shown in Figure 3. As observed ear- lier, random projection-based approaches achieve poor performance. When the features are not too noisy (as in rcv1 2 and rcv1 4), PCA-based methods are better than hdsl at compressing the space into very few dimensions, but hdsl eventually catches up. On the other hand, PCA suffers heavily from the presence of noise (dexter and dorothea), while hdsl is able to quickly improve upon the standard similarity in the original space. Finally, on all datasets, we observe that hdsl converges to a stationary dimension without overfitting, unlike pca+oasis which exhibits signs of overfitting on dexter and rcv1 4 especially.

CONCLUSION

In this work, we proposed an efficient approach to learn similarity functions from high-dimensional sparse data. This is achieved by forming the similarity as a combination of simple sparse basis elements that operate on only two features and the use of an (approximate) Frank-Wolfe algorithm. Experiments on real-world datasets confirmed the robustness of the approach to noisy features and its usefulness for classification and dimensionality reduction. Together with the extreme sparsity of the learned similarity, this makes our approach potentially useful in a variety of other contexts, from data exploration to clustering and ranking, and more generally as a way to preprocess the data before applying any learning algorithm.

}

  and use the resulting basis for the forward direction. This can be done in only O(M D) time and O(M D) memory and gives good performance in practice, as we shall see in the next section.

Figure 1 :

 1 Figure 1: Number of active features learned by hdsl as a function of the iteration number.

Figure 2 :

 2 Figure 2: Sparsity structure of the matrix M learned by hdsl. Positive and negative entries are shown in blue and red, respectively (best seen in color).

Figure 3 :

 3 Figure 3: k-NN test error as a function of the dimensionality of the space (in log scale). Best seen in color.

Table 1 :

 1 Complexity of iteration k (ignoring logarithmic factors) for different variants of the algorithm.

	Variant	Time	Memory
	Exact	Õ(T D 2

Table 2 :

 2 Datasets used in the experiments.

	Datasets Dimension Sparsity Training size Validation size Test size
	dexter	20,000	0.48%	300	300	2,000
	dorothea	100,000	0.91%	800	350	800
	rcv1 2	47,236	0.16%	12,145	4,048	4,049
	rcv1 4	29,992	0.26%	3,850	2,888	2,887

Table 3 :

 3 k-NN test error (%) of the similarities learned with each method. The number of features used by each similarity (when smaller than d) is given in brackets. Best accuracy on each dataset is shown in bold.

	Datasets identity rp+oasis pca+oasis diag-2	diag-1	hdsl
	dexter	20.1	24.0 [1000]	9.3 [50]	8.4	8.4 [773]	6.5 [183]
	dorothea	9.3	11.4 [150]	9.9 [800]	6.8	6.6 [860]	6.5 [731]
	rcv1 2	6.9	7.0 [2000]	4.5 [1500]	3.5	3.7 [5289] 3.4 [2126]
	rcv1 4	11.2	10.6 [1000]	6.1 [800]	6.2	7.2 [3878] 5.7 [1888]

Table 4 :

 4 Test error (%) of several SVM variants compared to hdsl. As in Table3, the number of features is given in brackets and best accuracies are shown in bold.

					Datasets svm-poly-2 svm-poly-3 svm-linear-2 svm-linear-1	hdsl
					dexter	9.4	9.2	8.9	8.9 [281]	6.5 [183]
					dorothea	7	6.6	8.1	6.6 [366]	6.5 [731]
					rcv1 2	3.4	3.3	3.5	4.0 [1915]	3.4 [2126]
					rcv1 4	5.7	5.7	5.1	5.7 [2770]	5.7 [1888]
		250
	Number of Selected Features	50 100 150 200
		0 0	2000	4000	6000	8000	10000
				iteration

In principle, any other convex and continuously differentiable loss function can be used in our framework, such as the squared hinge loss, logistic loss or exponential loss.

http://www.csie.ntu.edu.tw/∼cjlin/libsvmtools/datasets/

Note that the number of PCA dimensions is at most the number of training examples. Therefore, for dexter and dorothea, r is at most 300 and 800 respectively.

Note that hdsl uses roughly the same number of features as svm-linear-1 (Table4), but drawing any conclusion is harder because the objective and training data for each method are different. Moreover, 1-vs-1 SVM combines several binary models to deal with the multiclass

setting.5 Again, we were not able to run OASIS beyond a certain dimension due to computational complexity.

Acknowledgements

Most of the work was conducted when A. Bellet was affiliated with the U. of Southern California. This work was in part supported by the Intelligence Advanced Research Projects Activity (IARPA) via Department of Defense U.S. Army Research Laboratory (DoD / ARL) contract number W911NF-12-C-0012, a NSF IIS-1065243, an Alfred. P. Sloan Research Fellowship, DARPA award D11AP00278, and an ARO YIP Award (W911NF-12-1-0241). The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of IARPA, DoD/ARL, or the U.S. Government.


Isabelle Westeel

Cécile Malleret

Françoise Tayara

Aurélie Halipré

teaching and research institutions in France or abroad, or from public or private research centers.

The Open Access Week France event

Providing user guidance

A dedicated area for personalized workshops

With a customized Website and specific statistical tools, the Hal collections enhance the value of research laboratories, conference papers and research projects. The archives of Egyptologist Jean Vercoutter are the testimony of a researcher's life and were donated in 2012 to the Lille 3 University, where he taught Egyptology from 1960 to 1976.

First inventory of the

Special collections for research in Polish studies at Lille

  2011-2013: 10.000 items were catalogued and valued for researchers. Conservation and digitization projects are in progress for Polish newspapers published in France since WW1. Since 2014, the goal of the collaborative project Traditions in Antiquity in Lille and the Euroregion is to highlight and transmit national heritage by focusing on research, mediation and learning. The Hal ID is part of a researcher's digital identity. Being linked to the researchersGiven in 1911 by Auguste Angellier (1848-1911), the Laurenge-Angellier Donation (18.000 items) is part of the Heritage Collections of the University Library. He was the rst English language and literature teacher at the Arts Faculty of Lille, before being its Dean from 1897 to 1900. Its core areas are English history and studies, and most particularly the Irish 19th century literature and studies (completed by a collection of 4.000 items from the Patrick Rafroidi Library).

  ESPACE VIE ÉTUDIANTE BIBLIOTHÈQUE UNIVERSITAIRE CENTRALE -REZ DE CHAUSSÉEEntrée libre.Renseignements : hal-lille3.support@univ-lille3.fr La BU vous aide:-déposer dans HAL, -créer son identifiant, -générer une liste de publications. Venez avec votre ordinateur, on vous o re le café !


domain: phys spi chim

text:
Amanda J Neukirch

Wanyi Nie

Jean-Christophe Blancon

Kannatassen Appavoo

Hsinhan Tsai

Matthew Y Sfeir

Claudine Katan

Laurent Pedesseau

Jacky Even

Jared J Crochet

Gautam Gupta

Aditya D Mohite

Sergei Tretiak

Polaron stabilization by cooperative lattice distortion and cation rotations in hybrid perovskite materials

Keywords: organic-inorganic perovskite, polaron, photovoltaic, cation rotations

Hybrid organometallic perovskite (HOP) based solar cells have recently emerged as one of the most promising thin-film photovoltaic materials. [1][2][3] For example, solar cells prepared using low-temperature solution-based methods with methyl ammonium (MA=CH 3 NH 3 ) lead iodide, MAPbI 3 , have achieved an increase in efficiency from 3.5% to over 20% in about 5 years of effort. [4][5][6] The major draw of HOP materials comes from a combination of their high efficiency, low processing cost, and composition from earth-abundant and readily available elements. The unprecedented progress of HOP cell efficiency is indebted to research accomplishments in our understanding of photophysics of other electronically functional nanostructured materials such as organic semiconductors, metalorganic complexes, quantum dots, etc.

Despite the rapid development of organic-inorganic perovskite solar cells, there remain many questions about the fundamental photophysical processes taking place in these devices. For instance, HOP materials are known to have very low effective masses for electrons and holes, at ~0.15m 0 and 0.1m 0 respectively (here m 0 is the mass of a free electron) , [7][8][9][10] and very low electron-hole recombination rates, 11 but thin films exhibit modest charge carrier mobilities as compared to conventional semiconductors. 12,13 In order to capitalize on perovskite materials innate advantageous properties, a comprehensive understanding of the occurring phenomena is required. Just as HOP materials combine the efficiency of conventional semiconductor solar cells with the affordability of organic devices, lessons learned from both inorganic and organic materials are going to have to be integrated into one narrative.

Formation of polaron quasiparticles could prove to be very important for describing the charge carrier behavior in ionic and highly polar crystals such as perovskites. Lately many authors have proposed that the motion of the organic cations, specifically the reorientation of the dipole associated with CH 3 NH 3 + , can stabilize and localize charge carriers in the HOP lattice across many unit cells. [14][15][16][17][18][19] Electrons self-trapped over several lattice sites in a polarizable medium can be modeled in the effective mass approximation and are known as large polarons. [20][21][22] Large polarons typically appear in ionic or highly polar crystals (such as many III-V, II-VI semiconductors, alkali halides, and others), [22][23][24] where the Coulomb interaction between a conduction electron and the lattice ions results in an enhanced electron-phonon coupling. These quasiparticles have effective sizes much larger than the bond-length and travel through a lattice as free electrons but with an enhanced effective mass, which can explain the observed low electron-hole recombination rates 15,25 and modest mobility in perovskites 18 as compared to the values expected from the electron and hole effective mass. Polar fluctuations in HOPs can also support charge carrier separation which would contribute to long carrier lifetimes. 26,27 Additionally, large polarons have also been invoked to explain light induced changes to the HOP low-frequency dielectric response. 28 In contrast to the large polaron picture, when the binding energy is large, the effective mass approximation no longer works and the carriers are described as small polarons, which are spatially localized on a few lattice constants. 29 Small polarons play a very important role in organic semiconductors, [30][31][32] and transition metal oxides 24,33 by controlling charge carrier transport. They may also form in the HOP materials as indirectly evidenced by recent experimental observations of solar cell photodegradation dynamics. 34 In this case, small polarons would not contribute to the overall carrier mobilities, but are expected to create regions of charge accumulation which in turn cause performance degradation. Other than size, the most important difference between large and small polarons is in their charge transport: large polarons tend to have band-like transport, while small polarons usually undergo incoherent hopping transport.

Band and hopping transport are described by different theories and, therefore, in order to properly model transport in perovskite materials under different operating conditions, detailed studies on both small and large polaronic effects are needed.

Here we report the results from an in-depth hybrid Density Functional Theory (DFT) study of small polaron formation in HOPs. First, we invoke symmetry analysis in order to determine the lattices degrees of freedom coupled to the electronic system. Periodic boundary condition calculations were then performed to demonstrate how the energy levels and charge density of frontier orbitals were changed by polaron formation. Calculations of isolated clusters were used to explore charged state properties, localization of spin density and to quantitatively evaluate the upper bounds of reorganization energies, which altogether bespeak characteristic small polaron features. To commensurate our computational findings, we finally summarize recent experiments that provide indirect evidence for small polaron formation in a HOP solar devices. We find that upon full vibrational relaxation, small polarons are formed with sizes on the order of a lattice spacing, and binding energy deep into the band gap.

Symmetry considerations:

The building blocks of the hybrid lead halide perovskite structure are the Pb halide octahedra, PbX 6 4-(X=Halide), that are corner-shared to form a 3D crystalline lattice with a PbX 3 -stoichiometry. The cavities in the inorganic lattice are charge-balanced and filled with organic cations, such as MA=CH 3 NH 3 + . Thus, the crystalline structure of HOPs can be viewed as two alternating sublattices. The VBM (valence band maximum) is mainly the antibonding component of the hybridization between Pb s states and I p states, while the conduction band minimum (CBM) is dominated by the Pb p orbitals. Therefore, charge transport in HOPs occurs in the inorganic PbX 3-sublattice and the organic CH 3 NH 3 + sublattice acts as a medium that can modify the electrostatic landscape the charge carriers experience, leading to charge screening and localization. 18,35 Electron-phonon coupling is first investigated in the weak coupling regime assuming independent particles allowing group theory analyses. For the reference Pm-3m cubic space group of HOP, the possible intravalley non-polar deformation couplings of the VBM and CBM states only involve acoustic phonons, via a first order acoustic deformation potential (ADP) related to the strain tensor components. 36,37 Moreover, due to spin-orbit coupling in the CB, only local volumetric strain is expected to couple to VBM and CBM states and influence the first order ADP. This short range interaction is important and may yield in the strong coupling regime the so-called "acoustic polarons", which can trigger self-trapped electronic states. 24,38 On the other hand, none of the three lattice optical phonons can interact with the VBM and the CBM states via a conventional zero order non-polar optical deformation potential (ODP). 36,37 Before examining in detail the additional coupling between the VBM and CBM states and the molecular degrees of freedom, let's consider polar coupling mechanisms for phonons. The polar Frölich long range interaction between the delocalized CBM and VBM states and polar optical phonons (FOP) is an important carrier scattering process in hybrid perovskites as a result of dielectric increment between the high and medium frequency ranges related to longitudinal optical modes.

In the strong coupling regime, it yields the so-called "Fröhlich polarons". The classical piezoelectric electron-acoustic phonon polar coupling (PZA) vanishes in Pm-3m crystal structures. Simultaneous vanishing ODP and PZA electron phonon coupling mechanisms is a unique feature, which shows that halide perovskites define a specific new class of semiconductors.

In the parent cubic phase, the MA ion lies at the center of the regular cube formed by eight Pb 2+ ions. This requires the MA ion to be disordered in a certain way. 39 Discrete pseudospins (PS) models A, B, and C, Figure 1a-c, correspond to the thermally activated molecular tumbling occurring between the 6, 8, and 12 equivalent orientations from the center toward the faces ([100]), the edges ([110]), and the corners ([111]) of the cubic cell respectively. PS are simplified pictures of a more general approach involving continuous rotator functions. 40,41 Pseudospins and rotator functions allow describing molecular orientational fluctuations via non polar elastic and polar electric multipoles. Elastic multipoles span the same irreducible representations (IR) as the strain components, yielding a linear coupling with the lattice strain. 40,41 In the strong coupling regime, the formation of polaronic states is thus expected through deformation potentials and short-range interactions, between local volumetric strain, and local molecular arrangements once translational symmetry is lost, Figure 1 d,e. Static orientations of the cations correspond to a superposition of electric and elastic contributions. 41 Between knowledge of the charge density of band edge states and an understanding of symmetry allowed molecular arrangements lay the ground work for our understanding of the atomistic formation of polarons. Periodic boundary condition (PBC) DFT calculations were initially used to characterize the influence of local lattice deformations and eventual polaron formation. Early PBC calculations using DFT at the Generalized Gradient Approximation (GGA) level showed good agreement between theoretically determined band gaps and experimental measurements in the MAPbI 3 material. However, this agreement was due to the cancellation of errors between absence of spin orbit coupling (SOC) and many-body effects. 36,42 When both phenomena are included (e.g., with computationally expensive GW and SOC frameworks), the good agreement returns. 7,8 Hybrid functionals that include a fraction of the long-range exact exchange, such as HSE, 43 generally increase the band gap compared to standard GGA methods, thus improving agreement with experiment. Thus hybrid DFT has been proved to be numerically efficient method resulting in a reasonable accuracy when modeling perovskite materials. 8 In this work, for calculations involving PBC, we use the HSE and PBE0 functionals including SOC effects as implemented in the Vienna Ab initio Simulation Package (VASP) 44,45 (see Computational Details section in the SI).

Our modeling starts from the bulk structure for the low temperature orthorhombic phase of MAPbI 3 and optimizing until forces are below 0.02 eV/Å. 46 Based on the above mentioned symmetry analysis, suitable models can be constructed to illustrate the formation of polarons. To this end, to illustrate modification of frontier orbitals due to organic cation motions we consider a neutrally charged super cell of 192 atoms. Consistently with electrostatic interactions, the cationic MA's are rotated towards an electron residing at one of the Pb atoms mimicking the creation of a 'heavy' charge carrier or small polaron. The nearest neighbor MA cations are first manually rotated so that the dipoles are facing towards the chosen Pb atom. The system is then allowed to optimize for several optimization steps to reduce strain while retaining the direction of the MA dipoles.

Figure 2 shows the change in VBM and CBM states from the unperturbed system (a,b) and when the nearest 8 MA dipoles are rotated towards the Pb atom (c,d). In this case the charge density of the VBM remains highly delocalized away from the central Pb atom, while that of the CBM becomes localized on that Pb cation. This is indicative of the creation of a polaron for an excited electron in the CBM. The observed separation of charges is going to reduce radiative recombination between mobile and localized carriers. It should be noted that the charge densities remain unchanged with respect to cut-off energy, inclusion of SOC, van der Waal's interactions, and between functionals. In order to obtain a first estimate of the polaron stabilization energy, we calculated the change in band gap between the unperturbed system and the system where MA's have been rotated. Using both the HSE and PBE0 functional we got a band gap stabilization of about 160 meV in both cases. Isolated structure calculations: PBC, however, do not allow treatment of charged states due to spurious interaction of the charged cell with its replicas. We subsequently conducted simulations of finite structures similar to what has been done in the case of organic semiconductors. [30][31][32] When simulating charged systems, long range interactions are expected to play a significant role, requiring range-corrected functionals. For example, the hybrid exchange-correlation CAM-B3LYP functional (here fraction of the exact exchange varies from 20% to 65%) 47 in the presence of a polarizable solvent, reproduces correctly the polaronic effects in the organic semiconductors. [30][31][32] In this work, all calculations of HOP isolated clusters were performed using the CAM-B3LYP functional combined with the LANL2dz (for Pb and I) and 6-31G* (for N, C, and H) basis sets using Gaussian 09 software package 48 (see Computational Details section in the SI). Additionally, we include a polar solvent (ε = 78) via the conductor-like polarizable continuum model. 49 For comparison between the periodic system and isolated clusters, the calculations were also performed using the HSE model and these results are shown in the SI. There were several considerations when designing the isolated HOP structures. We started with lattice parameters that were derived from variable temperature powder x-ray diffraction experiments, taking the low temperature structure of MAPbI 3 as refined in a cubic lattice. 46 In order to minimize dangling bonds, the cluster should be terminated on all sides with the MAI (CsI). There must also be charge balance. Eqs. 1 and 2 ensure that the system is overall neutral

6! !" = 2! ! ! + ! ! ! (1) 2! !" + ! !" = ! ! ! + ! ! ! (2)
Here n Pb , n MA , and n I are the number of Pb, MA, and I atoms respectively. I 1 corresponds to an atom on the surface that is only bonded to one Pb atom, while I 2 corresponds to an I atom that is bonded to two Pb atoms. The cluster must be large enough to include at least 8 MA cations surrounding the central Pb atom. To fulfill all conditions, we end up with a cluster having the Performing calculations with both Cs and MA cations allows for separation between the effects of dipole rotation and volumetric strain. We will start by summarizing what takes place in the Cs-based system. The procedure for determining the polaron binding energy is to optimize the unperturbed neutral system. Once the geometry has been completely optimized, we add/remove an electron and perform a single point calculation to get the energy of the charged system in the neutral ground state geometry. We then compare this energy to a Cs-based system that has been geometrically optimized while charged. In order to quantify the polaron binding energy we calculate the two relaxation energies, the sum of which is the reorganization energy as, for example, enters in a standard Marcus theory

Nano Lett., 2016, 16 (6), pp 3809-3816; DOI: 10.1021/acs.nanolett.6b01218

describing electron transfer, [START_REF] Bao | Organic Field-Effect Transistors[END_REF] as shown in Figure S3 in the SI. The internal reorganization energy essentially quantifies variation of electronic energy due to the geometry changes when an electron is added or removed from the cluster. [START_REF] Marcus | [END_REF] The intramolecular reorganization energy is the sum of two components (λ 1 and λ 2 ) defined as

! !"#!$ = ! ! + ! ! = ! * -! + (! ! * -! ! ) (3)
Where λ 1 represents the change in energy between the neutral species in the optimized geometry of charged state (E * ) and the relaxed neutral species (E). The respective energy difference for the charged system at the neutral species geometry (! ! * ) and the relaxed charge geometry (! ! ), is represented by λ 2. The polaron binding energy is given by λ 2 . Using the CAM-B3LYP functional, we found that electrons have a binding energy of about 900 meV, which is much larger than the 270 meV binding energy for holes, Table 1. The results for the HSE functional are given in Table S1 and Figure S4 in the SI. As expected, binding energies for both holes and electrons decreased to 36 meV and 398 meV, respectively, owing to the lack of appropriate electronic exchange. Notably, semilocal models (e.g., GGA) and hybrid functionals with a small fraction of exact exchange (e.g., B3LYP, PBE0, HSE) are unable to describe polarons in other systems such as organic semiconductors. [30][31][32] As a reminder, all these numbers serve as an upper bound in polaron binding energy as the introduction of a finite system introduces quantum confinement effects. Moving on to the MAPbI 3 system, the procedure for finding the binding energy changes slightly. Once the neutral geometry has been optimized and a single point energy has been calculated for the charged versions, this energy is compared to a system where the MA molecules have been rotated (towards the central Pb atom for electrons and away for holes) and geometrically optimized in their charged state, as shown Figure 4 for the CAM-B3LYP functional. This ensures better starting structures for geometry relaxation of the charged states based on physical intuition. In the neutral geometry the hole electronic density is not in the center. This occurs because the geometry optimization spontaneously breaks the cluster symmetry and/or because the hole prefers to be centered on a halide which is not going to be centrally located. However in both cases the electron and the hole are significantly localized once the system is optimized in the charged state.

Electrons form a polaron binding energy of 1377 meV as given by the CAM-B3LYP model, Table 1. The magnitude of the polaron binding energy is likely over-estimated due to confinement effects, but it is clear that tighter polarons are formed when MA atoms are present.

The binding energy for holes is 578 meV. It should be made clear that not all of the binding energy comes from contributions of MA rotations as the system is also able to undergo volumetric strain during the optimization process. When these calculations are performed using the HSE functional (Table S1) the binding energies for holes and electrons become 33 and 880 meV respectively. The fact that both volumetric (Cs) and rotational strain (MA) can cooperate to small polaron formation indicates that the best way to reduce polaron formation, would be to use a larger organic cation with a smaller dipole moment such as Formamidinium (FA= HC(NH 2 ) 2 ).

These theoretical and experimental insights could explain why the current top perovskite solar cells are made with alloys containing FA. Experimental Evidence: We recently demonstrated the reduction in photocurrent, and subsequently device efficiency, in crystalline large grain size MAPbI 3 perovskite solar cells under constant solar illumination. The efficiency rapidly recovers to its original value after less than 1 minute in the dark. 34 Localized charged states strongly coupled to local structural lattice distortions and MA quasistatic configurations creating small polaron states, offer a nice explanation for the experimental observation of very slow degradation and fast self-healing on perovskite cell devices. In this picture, small polarons produced by light irradiation seed the formation of macroscopic charged domains preventing efficient charge extraction.

Conclusion:

Polaron quasiparticles frequently determine the charge carrier behavior in ionic and highly polar crystals such as perovskites. In this manuscript our computational study suggest small polaron formation in HOPs. A thorough investigation shows that localized charged states strongly couple to local structural lattice distortions and MA quasistatic configurations. Through symmetry arguments and DFT calculations using hybrid functionals, it was shown that MA dipoles can align with a charge centered on a Pb atom. Periodic boundary condition calculations demonstrate that frontier orbitals become localized after MA reorientation and that the band gap shows a stabilization in the 100's of meV. Moving on to charged isolated clusters, we found that electrons have a polaron binding energy much larger than holes, and are therefore able to predict that electrons are more readily trapped. Using the CAM-B3LYP functional we found that volumetric strain in a CsPbI 3 cluster creates a polaron with binding energy of around 900 meV.

In the MAPbI 3 cluster there is both volumetric strain and MA reorientation so binding energies are deeper at around 1300 meV. These computational insights allow rationalizing experimental studies showing that prolonged exposure to light leads to the accumulation of space charge, and yet there is no physical degradation or change in stoichiometry of the sample. Additionally, there is an evidence of a slowdown in MA rotation and formation of light-activated metastable trap states that are localized more than 500 meV deep within the band gap, 34 in excellent agreement with our numerical results. Such comprehensive computational studies are necessary in order to fully understand the complex behavior of perovskite materials and maximize electronic properties of these interesting materials to their full potential. For example, the fact that both lattice volumetric strain and rotational degrees of freedom can cooperate to create polarons, justifies superior performance of FAPbI 3 and Cs containing materials in solar cell applications. 6,55

Figure 1 .

 1 Figure 1. a-c Schematic representation of the totally symmetric IR of ! ! ! for the three pseudospins corresponding to [100], [111], and [110] respectively.39 Arrows indicate the orientations of C-N axis. d, e Artist view of a quasiparticle hole dressed by the interaction with neighboring rotating cations and travelling across the crystal.

Figure 2

 2 Figure 2 Localization due to rotation of MA cations. Crystal structures and charge densities including SOC of VBM and CBM a, b optimized experimental structure, c, d when 8 MA's are

  following stoichiometry: MA 54 Pb 27 I 108 (I 1 =I 2 =54). A cube fully terminated by MA (Cs) and I would require 64 MA molecules. In order to achieve charge balance between the ions and maximize symmetry the 8 corner MA's were removed as were two opposite MA molecules on opposite edges. FigureS1in the SI shows the resulting structures for both MA 54 Pb 27 I 108 and Cs 54 Pb 27 I 108 isolated clusters before and after geometry optimization. Despite the change in stoichiometry because the clusters originate from the bulk 3D geometry we will still refer to them as MAPbI 3 and CsPbI 3 .

Figure 3

 3 Figure 3 illustrates the spin density of electrons and holes in the optimized ground state and

Figure 3

 3 Figure 3 Hole (left) and electron (right) spin density distribution of optimized neutral (top) and charged (bottom) CsPbI 3 clusters calculated using the CAM-B3LYP functional.

Figure 4

 4 Figure 4 Hole (left) and electron (right) spin density distribution of optimized neutral (top) and charged (bottom) MAPbI 3 clusters calculated using the CAM-B3LYP functional.

Table 1

 1 Reorganization energies, in meV, for charged CsPbI 3 and MAPbI 3 clusters. Schematic representation of internal reorganization energy given in Fig.S3in SI.

	CAM-B3LYP		Electron			Hole
		λ 1	λ 2	λ	λ 1	λ 2	λ
	Cs	724	918	1642	365	274	639
	MA	510	1377	1887	528	578	1106

  6,[52][53][54] Recent reports on perovskite cells demonstrate that mixed alloys involving MA, FA, and Cs open the road toward further improvements, both for cell efficiency and stability.6,[52][53][54][55]

Nano Lett., 2016, 16(6), pp 3809-3816; DOI: 10.1021/acs.nanolett.6b01218

(Action Incitative, Défis Scientifique Emergents 2015). J.E. work was supported by the Fondation d'entreprises banque Populaire de l'Ouest under Grant PEROPHOT 2015. ACKNOWLEDGMENT AJN would like to thank Mikael Kepenekian for useful insight and discussion on the set up of the calculations.

Funding Sources

The work at Los Alamos National Laboratory (LANL) was supported by the LANL LDRD program (A.J.N., A.D.M, G.G and S.T.). This work was done in part at Center for Nonlinear Studies (CNLS) and the Center for Integrated Nanotechnologies CINT), a U.S. Department of Energy and Office of Basic Energy Sciences user facility, at LANL. This research used resources provided by the LANL Institutional Computing Program. LANL is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under contract DE-AC52-06NA25396. The work in France was supported by Cellule Energie du CNRS (SOLHYBTRANS Project) and University of Rennes 1

ASSOCIATED CONTENT

Supporting Information Available: Computational details and structures for Cs and MA isolated clusters before and after optimization, structures of isolated clusters, HSE charge densities for CsPbI and MAPbI 3 , schematic representation of reorganization energy, and experimental evidence for polarons. This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: ajneukirch@lanl.gov; serg@lanl.gov

Author Contributions

All authors have given approval to the final version of the manuscript.


Isabelle Mauger Perez

Olivier Cian

Isabelle Martin

French system for ETDs (and more) in a cartoon strip