Archetypal Analysis++: Rethinking the Initialization Strategy
Authors: Sebastian Mair
Jens Sjölund
Venue: Accepted by TMLR ()

Abstract: Archetypal analysis is a matrix factorization method with convexity constraints. Due to local minima, a good initialization is essential, but frequently used initialization methods yield either sub-optimal starting points or are prone to get stuck in poor local minima. In this paper, we propose archetypal analysis++ (AA++), a probabilistic initialization strategy for archetypal analysis that sequentially samples points based on their influence on the objective function, similar to $k$-means++. In fact, we argue that $k$-means++ already approximates the proposed initialization method. Furthermore, we suggest to adapt an efficient Monte Carlo approximation of $k$-means++ to AA++. In an extensive empirical evaluation of 15 real-world data sets of varying sizes and dimensionalities and considering two pre-processing strategies, we show that AA++ almost always outperforms all baselines, including the most frequently used ones.

Published in Transactions on Machine Learning Research (04/2024)
Archetypal Analysis++:: Rethinking the Initialization Strategy
Sebastian Mair
sebostian.mair wu.se
Uppsala University, Sweden
Jens Sjölund
jens. mnn.sjjounddki..
Uppsala University, Sweden
Reviewed on OpenReview: htips: // openreview. net/ met/ forum? id- KVUt 1M60RM
Abstract
Archetypal analysis a matrix factorization method with convexity constraints. Due to local minima, a good initialization is essential, but frequently used initialization methods yield either sub-optimal starting points or are prone get stuck in poor local minima.
In this paper, we propose archetypal analysis++ (AA++), probabilistie initialization strategy for archetypal analysis that sequentially samples points based their influence the objective function, similar to k-means++. In fact, we argue that k-means already approximates the proposed initialization method. Furthermore, we suggest to adapt an efficient Monte Carlo approximation of k-means++ to AA++. In an extensive empirical evaluation real-world data sets of varying sizes and dimensionalities and considering two pro- pro-processing strategics, we show that AA++ almost always outperforms all basclincs, including the most frequently used ones.
1 Introduction
Archetypal analysis (AA) (Cutler & Breiman, 1994] matrix factorization method with convexity constraints. The idea is to represent every data point convex combination of points, called archetypes, located on the boundary of the data set. Thus, archetypes can be seen as well-separated observations that summarize the most relevant extremes of the data. The convexity constraints also give archetypal analysis a natural interpretation. An example shown in Figure] 1
Archetypal analysis has been applied, e.g., for single cells gene expression (Thogersen al. 2013 Korem al.] 2015), bioinformatics (Hart et al.] 2015), apparel design (Vinéé et al.) 2015), chemical spaces of small organic molecules [Keller et al.] 2021), geophysical data (Black et al. 2022), large-scale climate drivers (Hannachi &
Convex Hull of X
Convex Hull of Z
Data Points x;
Initial Archetypes Zini
Trendafiov 2017 Chapman et al. 2022 and population
8e
Learning Path Z
genetics (Gimbernat-Maye
al.
2022
Learned Archetypes Z
Error
To improve the computation of archetypal analysis, vari- optimization approaches (Bauckhage & Thurau 2009
Mørup & Hansen
2012 Chen et al. 2014 Bauckhage 2015 Abrol & Sharma 2020 and approximations
ct al.
MMair et. al. 2017 Damle & Sun, 2017 Mei et al.) 2018
Figure 1: Archetypal analysis in two dimensions with à = randomly initialized archetypes {Z 24} shown in orange. The archetypes after optimization depicted in blue.
Mair & Brefeld 2019 Han et al. 2022) have been pro- posed. However, the carliest point of attack for obtain- ing a good solution is the initialization of the archetypes. Surprisingly, this has barely been investigated.

---

Published in Transactions on Machine Learning Research (04/2024)
In the original paper; Cutler & Breiman (1994) use random initialization, f.e., choosing points uniformly
at random from the data set, which was adopted by many others, e.g., |Egster & Leisch (2011 1; Seeth &
Eugster (2015); Hinrich et al. (2016): Hannachi & Trendafilov (2017); Mair et al. (2017): 0 Mei et al. (2018);
Krohne et al. (2019) Olsen et al. (2022) to name a few. Furthermore, Cutler & BBriiman (1994) state that a
careful lizitiiization improves the convergence speed and that archetypes should not be initialized too close
to each other.
This idea serves as an argument for using the FurthestPirst approach (Gonzalez 1985 Hochbaum & Shmoys, 1985), yielding well-separated initialization, e.g., for k-means clustering (Lloyd 1982). Inspired by Fur- FirstFirst, Morup & Hansen (2010; 2012) propose a modification for archetypal analysis called FurthestSum, which focuses on boundary points rather than well-separated points. Here, boundary points refer to points on the boundary of the convex hull of the data. Since then, FurthestSum has established itself as one of the default. initialization strategies for archetypal analysis which is used, e.g., by Thegersen et al. (2013)
Hinrich et al. (2016): Mair & Brefeld (2019); Abrol & Sharma (2020); Bock et al. [2022] Black al. (0022
Chapman el al. (2022 Gimbernat-Mayol et al. (022).
Despite its popularity, FurthestSum has also been criticized. For example, Suleman 2017 states that FurthestSum is prone to selecting redundant archetypes, primarily when many archetypes are used. Redundant archetypes lie in the convex hull of the already selected archetypes and thus do not lower the overall error. In addition, Krohne et al.] (2019) and Olsen et al. (2022) report better results with a random uniform initialization than with FurthestSum. A possible explanation is that FurthestSun's early focus on sub-optimal boundary points risks trapping the optimization of archetypal analysis in poor local minima.
Contributions. In this paper, we (i) propose archetypal analysis++ (AA++), an initialization strategy inspired by k-means++. Our theoretical results show that AA++ decreases the objective function faster than a uniform initialization in expectation. Furthermore, we (ii) argue that the k-means++ initialization can be seen as an approximation to our proposed strategy and that Monte Carlo approximation to the k-mcans++ initialization can be adapted for AA++ for more efficient initialization. Most importantly, we (ili) empirically demonstrate that our proposed initialization for archetypal analysís almost always outper- forms all baselines on 15 real-world data sets. Lastly, (iv) our extensive empirical evaluation serves as the first comprehensive evaluation of different initialization strategies for AA which is of independent interest.
2 Preliminaries
Before introducing archetypal analysis, we briefly revisit k-means clustering since we will build upon similar ideas. We use the following compact notation: [n]] = 21, n nz) for an E N.
k-means Clustering. Let X== fx,,;_ C Rd be a data set of n points in d dimensions, 2= fz Zk be a set of & cluster centers, and d(x, 3)² = mingez ||x - a/l3 be the minimal squared Euclidean distance from a data point x to the closest conter in Z. The k-means clustering problem has the following objective:
6x(2)= = Σ d(x, = Σ
min - 9ll..
XEX qEZ
XEX
Often, the cluster centers Z of k-means are initialized using the k-means++ initialization procedure (Arthur & Vassilvitskii] |007]] OOstrvsky et al.] 2013), which is guaranteed to yield a solution that is Ollog k)- competitive with the optimal clustering. The idea as follows. The first center is chosen uniformly random. Then, the remaining k- cluster centers are chosen according to a probability distribution where the probability of choosing a point x is proportional lo the minimal squared distance to the already chosen cluster centers, i.e., p(x) x d(xx 2) The procedure is outlined in Algorithm 3 which can be found in Appendix B
Archetypal Analysis. Let = {x:] C Rd be a data set consisting of n E N d-dimensional data points arranged as rows in the design matrix X E Raxd The idea of archetypal analysis (AA) (Cutler & Breiman 1994) is to (approximately) represent every data point x; as convex combination of k E N archetypes
2

---

Published in Transactions on Machine Learning Research (04/2024)
Z {Z1. (21...., zz), i.e.,
x R a Z, a =1=1, a, 2
(1)
where the matrix E Rkd contains the archetypes as rows and the vector a, E R defines the weights for the ith data point. Here, 1 denotes the vector of ones and a, 2 0 meant element-wise. The archetypes zj (j [k:)) themsolves are also represented (exactly) as convex combinations, but of the data points, i.e.,
bJx, = bj b020,
where by E R" is the weight vector of the jth archetype. Let A E R and B R be the matrices consisting of the weights a, (i E [n]) and by (j E [k]). Then, archetypal analysis yields an approximate factorization of the design matrix X as follows
X ABX = AZ,
(2)
where BX R is the matrix of archetypes. Due to the convexity constraints, the weight matrices A and B are row-stochastic. The weight matrices A and B are typically determined by minimizing the approximation error in Frobenius norm, resulting in the following optimization problem
minimize A.B
X A
subject to
A1 = A20 and B1= B20
This can be equivalently expressed minimizing the sum of projections of the data points the archetype- induced convex hull as follows
X AZI
Σ
min ||x q -12
(4)
9Econv (2)
XEX
where conv(2) refers to the convex hull of the set 2. An example of archetypal analysis and the projection generalized low-rank problem (Udell et al., 2016), which is biconvex but not convex. Because it is biconvex, local optimum can be found via an alternating optimization scheme introduced Cutler & Breiman (1994), such as the standard one outlined in Algorithm dependent on the initialization.
errors is depicted in Figure [ The optimization problem
! which can be found in Appondix
B However, the quality of such local optimum is directly
Archetype Initializations. Popular ways initializing the archetypes is by using data points chosen uniformly at random (called Uniform) or the FurthestSum procedure (Morup & Hansen, 2010 2012), but we also introduce FurthestFirst. Originally proposed for the metric k-center problem, the FurthestFirsl algorithm [Gonzalez] 1985 Hochbaum & Shmoys 1985) selects the first center/archetype uniformly at ran- dom and iteratively adds the point that is furthest away from the closest already selected center center/archelle /archetype. Formally, the index of the next point is j juext arg max,e (mingez ||x; - -9ll2). Specifically for archetypal analysis, Merup & Hansen (2010 2012) propose a modification of Furthest First called FurthestSum, which
sums over the distances of the already selected points, i.e., juext arg max;e (Eqez ||x, - -xll)) To increase the performance, the first point, which was chosen uniformly at random, is usually discarded in the end and replaced by new point chosen via the criteria outlined above.
Although not designed as an initialization strategy, we also mention A.corrsee Mair & Brefeld 2019) as it was used for the initialization of archetypes by Black et al. (2022) and Chapman et al. (2022 A coreset is a small subset of the original data set that allows for a more efficient training of archetypal analysis. The selection probability of every data point is given by p(x) x ||x pll2, where μ is the mcan of the data set X.
3 Archetypal +
The idea of our proposed archetypal analysis++ (AA++) initialization procedure is very similar to the one of k-means+++. We begin by choosing the first archetype zi uniformly at random. The second archetype
3

---

Published in Transactions on Machine Learning Research (04/2024)
Algorithm Archetypal + Initialization
Input: Set of n data points X, number of archetypes k
Output: Initial set of k archetypes Z
Sample an index i unilormly al random from [n]] i.e., using p(i) = nn
4: Append X; to 2
5: while || < k do
6: Sample using p(t) x mingeconv(2) ||x; - -912
7: Append x; to
end while
k=1
k=2
Unifom
MSE=39.40
MSE=8.55
MSE=4.60
MSE=4.55
SurthenSSum
MSE-59.11
M5E-15.03
MSE-0.09
MSE=0.08
MSE=39.44 M5E=12.54
MSE=0.54
MSE=0.31
Figure 2: A comparison of Uniform, FurthestSum, and the proposed AA+ when consecutively initializing k = 4 archetypes. MSE denotes the mean square error, i.e., Equation multiplied by n
is chosen according a distribution that assigns probabilities proportional to the distances from the first archetype, i.e., p(x) x ||x 211.2- The remaining
archetypes are chosen according to probability distribution where the probability of choosing a point x is proportional to the minimum distance to the convex hull of the already chosen archetypes, i.e., p(x) x minn | qll.. This procedure is outlined in Algorithm
With every additional point sampled by AA++, the convex hull of the initialized factors, é.e., conv expands. This because selecting point outside the convex hull of Z, i.e., a point in (x E X x & conv(2)), is by dcfinition not contained the convex hull of and hence expands it. In contrast, a point within the convex hull of Z would not increase its volume but has zero probability of being selected by AA++ Note that uniform initialization also assigns probability mass to points in the interior of conv(Z). Selecting new point can make a previously selected point redundant, i.e., it then lies in the convex hull of previously selected archetypes and does not help to increase its convex hull. However, this is not a problem in practice since redundant. archetypes can be recovered during optimization.
Note that the optimization problem in line the archetypal analysis implementations can be re-used, simplifying the implementation AA++. Besides, AA++ has no hyperparameters, and line is trivially parallelizable
Algorithm 1 is the same as in Equation 1 Thus, parts of
4

---

Published in Transactions on Machine Learning Research (04/2024)
Figure _ depicts how Uniform, FurthestSum, and our AA++ strategy initialize 2 on synthetic data set. The gray projection errors can be seen as selection probabilities for AA+ In contrast, the selection probabilities of Uniform (not depicted) are not data-dependent,, and FurthestSun sums the distances to the already selected archetypes up, i.c., not the depicted distances. The sum of all depicted projections resembles Equation [ and we provide the Mean Squared Error (MSE), i.e., Equation (4) multiplied by n n-1, per iteration k.
In terms of MSE, FurthestSum works slightly better than +A++ on this synthetic data set. However, this is no longer the case with higher-dimension real-world data, our empirical evaluation later shows.
Theoretical Analysis. We first show that by adding a new archetype to the set of archetypes in the while loop of AA++ in Algorithm T the objective function is guaranteed to decrease.
Lemma 3.1. Let ||X - AZIl > 0, i.e., there are points yielding projection errors. Then, adding a point x X Z the set of archetypes Z according to AA++ (Algorithm 1 is guaranteed to decrease the objective function.
Note that a uniform initialization does not have this guarantce since it also assigns probability mass to points that do not yield a reduction. We now show that, in expectation, AA++ reduces the objective function more than the uniform initialization; equality is just met if all candidate points have equal projection distance.
Proposition 3.2. Let [Z] 2 1 and X - AZI > 0. Per iteration, AA++ (Algorithm i yields a reduction of at least as much as the reduction for the uniform initialization in expectation. The reduction is equal only if every point in X Z has the same projection onto the current set of archelypes 2. Furthermore, using uniform initialization cannot reduce the objective function more than AA++ in expectation.
The proofs of both statements can be found in Appendix
Complexity Analysis. The proposed initialization strategy outlined in Algorithm 1 selects the first archetype uniformly at random. The remaining k - archelypes are chosen according to a probability proportional to the squared distance between the candidate point and the convex hull of the already chosen archetypes. To compute this projection, quadratic program (QP) has to be solved. Thus, the complexity of Algorithm D is O(n- (k 1)- QP). It depends not only on the size of the data set and the number of archetypes & to be initialized but also on the complexity of solving the quadratic program O(PP), which is often cubic in the number of variables & (Goldfarb & Liu 1990).
4 Approximating the Archetypal + Initialization
The proposed inilialization procedure AA++ has to solve n (A: 1) quadratic programs, which can be time--oonsuming Thus, we propose two strategies to approximate the initialization procedure. The first one approximates the distance computation (circumventing the computation of quadratic programs), and the second one approximates the sampling procedure (reducing n).
Approximating the Distance Computation. Mair & Brefeld (2019) show that the objective function of k-means upper bounds the objective of archetypal analysis, which is due to the per-point projections, i.e.,
min qEconv(2)
||x - 1 < min |x -91
QE2
for set of clusters/archttypn Z. Since the computationally most expensive operation in the proposed Algorithm D is the projection onto the convex hull in line 6, it can be approximated by the distance to the closest point within the already chosen archetypes. See Figure 3 for an example. Note that the distance is then always over-estimated, and even points within the convex hull of the already chosen points might have a distance, although the projection should have a length of zero. This depicted in Figure 3 for the red point. Following this approach boils down to the k-means initialization procedure. Hence, the new complexity is Onn k- d), thus avoiding the cost of solving the QP.
5

---

Published in Transactions on Machine Learning Research (04/2024)
Convex Hull of
Approximated Distance Initialized Archetypes
True Distance
Figure 3: Approximation of the distance function in two dimensions. The true distance of the green point is depicted using solid line whereas the approximation is shown as a (larger) dashed line. The red point has no distance to the convex hull, but the approximation yields a positive distance.
Approximating the Sampling Procedure. Another approach approximate the sampling proce- dure by not considering all data points, which is especially beneficial in large-scale scenarios. To initialize k-means++ in sublinear time, Bachem el al. (2016) leverage Markov Chain Monte Carlo (MCMC) sam- pling procedure. This procedure is based on the Metropolis-Hasting algorithm (Hastings) 1970 with an independent and uniform proposal distribution. We adapt this idea for AA++ as follows. The first archetype is still sampled uniformly at random. For every following archetype, a Markov chain length m < n is constructed iteratively. We begin by sampling an initial point x,. Then, in every step of the chain, we sample candidate point x, and compute the acceptance probability, which is given by
m 2) |x, 1 Xx: q//
min min(1,
mill
With probability म, we update the current state from X to Xj, otherwise we keep X After m steps, we add the current x; as the next archetype to the set of initial archetypes 2. This approach is summarizod in Algorithm 2 and the complexity of this strategy O(m k QP).
Bachem al. (2016) also provide thcoretical result that bounds the error in terms of the total variation distance of the approximate sampling distribution to the true sampling distribution.
Algorithm AA++ Monte Carlo Initialization
1: Input: Set of n data points X, number of archetypes k, chain length m
2: Output: Initial set archetypes Z
Sample an index uniformly at random from [n], i.e., using p(i) n!
Append x, to Z
while |Z| k do
6: Sample i uniformly at random from [n], i.e., using p(i) n
7: Compute the distance to the convex hull, i.e., of ming |xx; 911
8: for 1 2.3...,, do
Sample j uniformly at random from [n], i.e., using p(j) = -1 Compute the distance the convex hull, i.e., df ming ||x, -/1. Sample r uniformly at random from (0, 1)
10:
11:
if d? = or d then
12:
13:
Set and df d3 14: end if
15: end for
16: Append x; to 2
17:
end while

---

Published in Transactions on Machine Learning Research (04/2024)
Here, |lp rllrv denotes the total variation distance between two distributions p and q which is defined
lp <||Tv = 2/2 Σ [((x) - -x))]=
XEQ
where is a finite sample space on which both distributions are defined on. The following bound on the error shows that the longer the chain length m, the smaller the error c.
Theorem 4.1 (Adapted from Bachem et al. (2016)). Let k > 0 and 0 < 1. Let p++ be the probability distribution over Z defined by using AA (Algorithm defined by using AC+ (Algorithm 2
7
and PACME be the probability distribution over 2
Then,
||MMMMC P++IITV Se
for a chain length m = O( log ) where
d(x, 2)² Ex'ex 2))
max TI 2CR XEA
max
and d(x, 2)² 3:0=mim ( ||x 411..
Note that y is property of the data set.
5 Experiments
Data. We use the following seven real-world data sets of varying sizes and dimensionalities. Additional eight real-world data sets, often smaller in size and dimension, are considered in Appendix D
The California Housing (Pace & Barry 1997] data set has n = 20, 640 examples in d = dimensions. A first large data Covertype [Blackard Dean 1999), which consists of n 581, 012 instances in d= 54 dimensions. FMA Defferrard et al.] 2017) a data set for music analysis that considers n 106, 574 songs reprosented with d 518 features. It thus serves as data set of higher dimensionality within evalua- tion. KDD-Proteir 2 has 145, 751 data points, each represented with d = 74 dimensions measuring the match between a protein and a native sequence. The data set Pose is a subset of the Human3. 6M data set (Catalin Tonescu 2011 Tonescu et al. 2014). Pose was used in the ECCV 2018 PoscTrack Challenge and deals with 3D human posc estimation Each of the n = 35, 832 poses is represented as 3D coordinates of 16 joints. Thus, the problem is 48-dimensional. The data set RNA (Uzilov et al., 2006) contains n = 488, 565 RNA in- put sequence pairs with d = 8 features. Another larger data set we use is subset of the Million Song Dataset MBartin-Mahieux et al. 2011). which is called Song. It has = 515, 345 data points in d = 90 dimensions.
Data Pre-processing. We apply pre processing to avoid numerical problems during learning and consider two different approaches: (i) CenterAnddoarSeale which centers the data and then divides the data matrix by its largest element, and (ii) Standardizalion which centers the data and then divides every dimension by its standard deviation. Results on CenterAndMoozSolle are presented in the main body of the paper, while results on Standordization are presented in Appendix G for completeness.
Baselines and Approximations. We compare our proposed initialization strategy AA++ against Uniform subsample of all data points, FurthestFirst (Gonzalez] 1985 Hochbaum & Shmoys, 1985). FurtheslSum (Morup & Hansen 2010 2012), and A Acoreset (Mair & Brefeld] 2019). In addition, we evaluate our two proposed approximations of AA++. The first. approximation, which is equivalent to k-means++, should thus not be seen as a baseline. Note that by k-means we only refer to the initialization strategy, i.e., we do not run k-means afterwards. The second approximation uses a Markov Chain and called AA MC.. We evaluate 1%, 5%, 10,, and 20% of the data set sizes as chain lengths m.
https://github... http:/tttuub.ooommee/ee//eee//// www/ndeff//ma.aam
http://osmot.ce wtw.ccccccc.l....oo. wwt/kiddcppdsttaeetet... html
http://vision inar. www.:rrr//aaa.eemm wwt.eaachhaaaaaaaaaeaaaaaaaeaaaa chaleeeeeeeaaleeeeoeeeoeoeeee .ptpeepppppppppppp

---

Published in Transactions on Machine Learning Research (04/2024)
Unfom
FudhssFros
Sur
Alcorsct
MMars:
M11
--- AAr ICC 20%
AA1 RC 10%
AA IMC 5%
A4 IMC 15
California X-15
Califoria C-25
Califorii R-50
California A-75
California k-100
107
10"
MSE
20-0
7
107"7
10"
Coveraype 1--15
Covertype k-25
Covertype 6-50
Covertype Aw75
Covertyoe k=100
111
N00E
10"
FMA 6-15
2-25
FMA ビ-50
FMA
FMA
11
11
KDD-P Protein 15
DDD-Protein 2-
KDD-Protein Protein E-50
DDD-Pratein 15
KDD-Protein 1=100
MSE
10
10
Pose 2-15
Poss 2--5
Pose k-50
Pose k-75
Poss 1o-300
11
4
1.5 100
8x 10-1
3x -0"1
MSE
10
10
2x 10
4x 11-1
3
2 10m1
3xx 11111
RNA A-15
RNA k-25
RNA k-50
RNJA A-75
RNA k-100
111
100
MSE
10
10
10"
Song km
Song k-25
Song im30
Song
Song k-100
10"3
10
10
10
1x 10
10
10
3xx 10
1rh.
10
15
20
25
14L
13
15
20
25
30
44L
10
15
20
25
112.
5
10
15
20
25
165.
10
15
20
25
30
Iarabone of AA
herat oma of AA
harattons of AA
Itaratione AA
Iberstinne of AA
Figure& 4: Results on California Housing, Covertype, FMA, KDD-Protein, Pose, RNA, and Song.
8

---

Published in Transactions on Machine Learning Research (04/2024)
Best Initialization
Best Overall
Median Initialization
Median Overall
1=15
2
0
11
3
1
9
1
10
2
9
14
12
2=
0
3
0
0
10
0
2
0
2
11
3
0
0
11
2
0
0
0
12
10
5=
0
3
0
1
0
11
1
0
0
0
13
0
2
1
0
0
12
2 0
0
0
12
8
6
75
1
0
0
13
1
0
1
1
11
1
1
0
0
13
0
0
0
14
100
0
14
1
12
1
13
1
12
Uniform
FurthestSun
AAcoreset
k-Means+ +-+
AA+-
Uniform
FurthestFirst
FurthestSum
AAcoreset
A-Means++
AA++
Uniform
Furthesst
First
+-+
AA+-
Uniform
FurthestSun
Accorsset +-+ AA+-+
k-Means-+
AA
Figure 5: Aggregated statistics over 15 data sets (seven data sets from above and eight data sets from the appendix). Each table shows how often each initialization method yields the best result for various choices of k under different settings. Best refers to the lowest error of a single seed and median refers to the median over many seeds. We report on the performance after initialization and overall during the optimization.
Setup. For various numbers of archetypes k E {15, 25, 75, 100), we initialize archetypal analysis according to each of the baseline strategies and compute the Mean Squared Error (MSE), í.e., the objective in Equation 4 normalized by n-', In addition, we perform a fixed number of 30 iterations of archetypal analysis based on those initializations. Here, we use the vanilla version of archetypal analysis according to
Cutler & Breiman (1994). We compute statistics over 30 sceds, except for larger data sets (n > 500, 000 d >> 500) for which we only compute 15 seeds and report on median performances and depict the 75% and 25% quantiles. Furthermore, we track the time it takes to initialize the archetypes as well as the time each iteration of archetypal analysis takes. Details on the implementation are provided in Appendix C All experiments run on an Intel Xcon machine with 28 cores with 2.60 GHz and 256 GB of memory.
Performance Results. Figure depicts the results on California Housing, Covertype, FMA, KDD- Protein, Pose, RNA, and Song. Although being the most commonly used initialization methods, Uniform (red line) and FurthestSum (blue line) often yield the worst initializations, which can be seen by a short straight line on the left-hand side. As expected, the error decreases during the first 30 iterations of archetypal analysis. However, there is frequently a significant performance gap between Uniform when compared to AA++ (green line), especially for California Housing, Covertype, FMA, and RNA. The same gap is visible for FurthestSum on all data sets. The distance approximation which equivalent to k-means++ (orange line) sometimes performs similarly to AA++ but is almost always better than Uniform. As for the Monte Carlo approxima- tions of AA++ using 1%, 5%, 10%, and 20% of the data set size as chain lengths m, we can see that they approximate AA++ sufficiently well on these data sets, except on California Housing and RNA, where the smallest chain length struggles. Overall, our proposed AA++ initialization performs best in almost all cases.
This is confirmed by the aggregated statisties depicted in Figure 5 initialization method (MCMC approximations excluded) performs best on all evaluated data sets. In total, we consider 15 data sets: all previously introduced data sets and eight data sets, which are discussed in the appendix. Thus, 15 is the highest and best number that can appear in each of those lables. Within Figure 5 best refers to the lowest error among all seeds, and median refers to the lowest median over all seeds per method. Besides, initialization only considers the performance after initialization, whercas overall reflects the best performance across all 30 iterations of AA, which is lypically at the last iteration. Once again, AA++ clearly yields the best results in all scenarios, especially for larger values of k. Note that we show k-means separately for completeness, but it should be seen as an approximation of AA++. For completeness, we also report on aggregated statistics on worst performances in Figure which can be found in the appendix.
Each table shows how often each
Timing Results. The incrcase in performance comes at a cost. Figure G depicts the time needed initialize the archetypes on the California Housing and Covertype data sets for various choices of k. As expected, the fastest-performing method is Uniform since it uses no information about the data except the
9

---

Published in Transactions on Machine Learning Research (04/2024)
Uniform
AAcoreset
1x AA Iter.
AA+ MC 10%
FurthestFirst
Means+
&x Iter.
AA+ 5%
Fur
++++
AA MC 20%
AA+
MC 1%
California
Covertype
Initialization Time Seconds
10
10a
10
10-1
10-2
10-3
25
100
15
50
75
100
h
k
Figure 6: The median time it takes to initialize archetypal analysis.
number of data points. AAcoreset slightly slower as needs to compute the mean of the data as well all distances the mean. FurthestSum is slower than Uniform and AAcoreset, and the proposed approach AA++ is consistently the slowest. However, note that initializing k points is still much cheaper than rumning & AA iterations (black line). Using the distance approximation which is equivalent to k-means- is much faster than AA++ and takes approximately the MCMC-based approximation also drastically reduces the initialization time of AA++. On Coverlype, AA++MC using 1% of the data points as the chain length m takes approximately as much time as FurthestSun, and other chain lengths m lower the initialization time to the time of k-means+
much time to initialize as FurthestFirst. Using
Influence of the Pre-processing Scheme. A surprising finding that FurthestSum is especially sensi- tive the propprocessig scheme. Standardization clearly often degrades the performance FurthestSum but does not affect + and its approximations as show in Appendix G
6 Discussion & Limitations
The idea of archetypal analysis++ is closely related to k-means++. Thus, it is natural to expect that a similarly strong thcoretical guarantec, i.c., E Ejok-means++] S O(log k)oOPT, can be derived for AA++ as for A-means++. Lloyds algorithm (Lloyd] 1982) consists of two rather simple steps that are iterated until convergence. First, within the assignment step, the distances between data points and cluster centers are computed, and each data point is assigned to its closest cluster center. Second, per cluster, the cluster centers are updated by computing the mean of all points assigned to this cluster. Note that those computations are relatively simple, and the second step can be performed in closed-form. Due to this simplicity, many theo- retical analyses on k-means (including k-means++)) exist. In contrast to that, archetypal analysis is a more challenging problem, not only computationally due to the optimization problems that do not admit a closed- form solution (ef. Algorithm 4), but also theoret to analyze. For k-means, Arthur Vassilvitskii (2007) exploit having closod-form solutions for the optimal cluster centers, and that only points in a cluster influence closed-form
the position of their cluster center. This is not possible for archetypal analysis since there constructing convex hull. Thus,
solutions, and there is an interdependence between archetypes since we believe similar strategy to derive theoretical guarantecs for AA+ is unfruitful. However, note that neither of the state-of-the-art initializations for archetypal analysis come with any theoretical guarantees.
Undoubtedly, data-indens indepeededddd initialization strategies like Uniform that chooses points from the data set uniformly at random are the fastest way initialize model. Any data-dependent strategy will be more costly in terms of initialization time. One drawback of the proposed AA+ initialization strategy is that it is computationally more involved than its baselines. More specifically, it has to solve quadratic program
10

---

Published in Transactions on Machine Learning Research (04/2024)
per archetype to be initialized except the first one. This renders the initialization of k archetypes typically to be more expensive in terms of time than a single iteration of archetypal analysis but faster than k iterations, cf. Figure 0 As a remedy, we proposed two approximations of AA++, which both decrease the initialization times to those of k-means++ or a single iteration of archetypal analysis. Thus, allowing for a practical usage of our proposed initialization stralegy.
We might expect that the total runtime (30 iterations of archetypal analysis) is mainly determined by the different initialization efforts, i.e., archetypal analysis is being the slowest when initialized using A However, this is not true. As seen in, e.g., Figure [0] the curves start indeed at different times, and AA- is usually the most computationally demcaning initialization approach. Nonetheloss, the total time of AA-+ is not always the slowest, see, e.g., Covertype in Figure 10 Note that the figure depicts the performance after initialization, followed by 30 interpolated points that correspond to the performances after 30 iterations of archetypal analysis. Per iteration, lines 5-8 in Algorithmmaae computed. Specifically, line involves solving n non-negative least-squares problems to obtain the mixture matrix A; line computes the archetypes Z given A and the data matrix X by solving a system of linear equations; line 7 involves solving & non-negative least- squares problems to obtain the mixture matrix B; and line 8 computes the archetypes Z given B and X by computing matrix-matrix product. There are several factors that determine the runtime of cach and every iteration. First, to solve the non-negative least-squares problems, we utilize NNLS [Lwsson & Hanson 1995). which is an active set method that, depending on the specific inputs, needs a different number of internal iterations, and thus, the runtime of every non-negative least-squares problem varies. Second, instead of solving a system of linear equations in line 6, we solve a least-squares problem via np. ap..linag..... lssss for cases in which A is not full-rank. Both factors can influence the runtime of an iteration of archetypal analysis.
Due to the bi-convex loss and the alternating optimization scheme, we expect the loss to decrease in every iteration. However, on rare occasions, e.g., RNA and Song in Figure E the loss slightly increases. We conjecture that those cases are related to numerical instabilities. For example, instead of solving a system we solve a least-squares problem via np. ap.linalg... 1stsq for cases
of linear equations in line 6 of Algorithm ! in which A is not full-rank. This might have an influence on the positioning of the archetypes Z.
7 Related Work
Other variants exist besides the classical archetypal analysis (Cutler & Breiman 1994b. Moving archetypes
Cutller & Stone 1997) define AA for moving targets. There are adaptations of archetypal analysis for missing data (Epilanio et al. 2020) and interval data (D'Esposito et al.] 2012). and probabilistic archetypal
analysis (Seth & Eugster 2016) rephrases the factorization problem in a probabilistic way. De Handschutter
ct al. (2019) relax the convexity constraints of AA. More recently, approaches basod on deep learning have : van Dijk et al. (2019); Keller et al. (2021). For those, the initialization of archetypes is irrelevant, as considered in this paper. Furthermore, archetypoid analysis ςVinué et al., (2155) restricts the archetypes to be data points instead of convex combinations of data points. Hence, the idea is similar to k-modoids (Kaufman & Roussccuw
been considered Keller et al. 2019
1990) and intends to add even more interpretability. Note that our proposed AA++ can also be used as an initialization for archetypoid analysis.
Suleman 2017) stresses that an improper initialization of archetypal analysis is problematic and that the FurthestSum method is prone to selecting redundant archetypes, especially when having many archetypes. Nascimento & Madaleno (2019) consider an anomalous pattern initialization algorithm for initializing archetypal analysis. However, their focus was on inferring the number of archetypes l and then using archetypal analysis for fuzzy clustering.
Less common initialization strategies for archetypal analysis include rumning k-means clustering, as used by Han et al. (2022), and utilizing a coresct (Mair & Brefeld 2019). as used by Black et al. (202)) and Chapman el al." (022). Note that the coresel was proposed as a way to condense the data set into a smaller set for a more efficient training of archetypal analysis rather than as a way for initializing it. However, we still include it as a bascline for completeness in our experimental evaluation. As for k-means, we drew inspiration from the k-means++ strategy but do not consider initialization strategies for k-means as suitable basclines for archety- pal analysis. We rather focus on initialization strategies that are frequently used for archetypal analysis.
11

---

Published in Transactions on Machine Learning Research (04/2024)
problem, several initialization techniques based on randomization, other low-rank decompositions, clusterings, heuristics, and even learned approaches (Sjölund & Bånkestad] 2022) are used. A summary of various NMF initialization methods is provided by Esposito (2021). Among the considered baselines in this paper, all are randomized. However, FurthestSumn can be seen as a heuristic algorithm, and Furthest First and k-means+++ can be classified as (ini- tializations of) clustering algorithms. The proposed AA++ approach randomized just as Uniform, with the difference that the selection probabilities in AA++ are data-dependent while Uniform data-independent.
For the related non-negative matrix factorization (NMF) Lee & Seung 1999)
8 Conclusion
We introduced archetypal analysis++ (AA++) as a new initialization method for archetypal analysis that is inspired by k-means+++ The proposed method does not have any hyperparameters and straightforward to implement,
allows the reuse of already existing subroutines of any implementation of archetypal analysis. Furthermore, proposed two approximations of AA++ that are faster to compute. The first one approximates the distance computation and is equivalent to the k-means++ initialization method. The second one is an MCMC approximation of the sampling procedure of AA++. which reduces the number of quadratic problems that have to be computed. We showed that both approximations work reasonably well and lower the initialization times to competitive levels. Empirically, we verified that, for two pre eve-processing schemes, AA++ almost always outperforms all basclines, including the most frequently used oncs, namely Uniform and SurthestSum, on 15 real-world data sets of varying sizes and dinensionalities.
Acknowledgments
This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. The authors thank the anonymous reviewers, Samuel G. Fadel, Ahcène Boubekki, Paul Häusner, and Zheng Zhao for valuable comments on earlier drafts.
References
Vinayak Abrol and Pulkit Sharma. A geometric approach to archetypal analysis via sparse projections. In Proceedings of International Conference on Machine Learning, pp. 42 51. PMLR, 2020.
David Arthur and Sergei Vassilvitskii. k-means-4: The advantages of careful seeding. In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA '07, pp. 1027-1035. Society for Industrial and Applied Mathematics, 2007.
Olivier Bachem, Mario Lucic, S Hamed Hassani, and Andreas Krause. Approximate k-means++ in sublinear time. In Thirtieth AAAI Conference on Artificial Intelligence, 2016.
Christian Bauckhage and Christian Thurau. Making archetypal analysis practical. In Joint Pattern Recog- nition Symposium, pp. 272-281. Springer, 2009.
Christian Bauckhage, Kristian Kersting, Florian Hoppe, and Christian Thurau. Archetypal analysis as autoencoder. In Workshop New Challenges in Neural Computation, volume 2015, pp. 8-16, 2015.
Connor Beck, Anja Kunze, and Dominique Zosso. Archetypal analysis for neuronal clique detection in low- rate calcium fluorescence imaging. In Proceedings of the 44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society, pp. 162-166. IEEE, 2022.
Thierry Bertin-Mahieux, Daniel P.W. Ellis, Brian Whitman, and Paul Lamere. The million song dataset. In Proceedings of the 12th International Conference on Music Information Retrieval, 2011.
Amanda S Black, Didier P Monselesan, James S Risbey, Bernadette M Sloyan, Christopher C Chapman, Abdelwaheb Hannachi, Doug Richardson, Dougal T Squire, Carly R Tozer, and Nikolay Trendafilov. Archetypal analysis of geophysical data illustrated by sea surface temperature. Artificial Intelligence for the Earth Systems, :3)=220077 2022.
12

---

Published in Transactions on Machine Learning Research (04/2024)
Jock A Blackard and Denis J Dean. Comparative accuracies of artificial neural networks and discrimi- nant analysis in predicting forest cover types from cartographic variables. Compulers and Electronics in Agriculture, 24(3):131- 151, 1999.
Thomas F Brooks, Stuart Pope, and Michael A Marcolini. Airfoil self-noisc and prodiction. 1989.
Cristian Sminchisescu Catalin Tonescu, Fuxin Li. Latent structured models for human pose estimation. In International Conference on Computer Vision, 2011.
Chih-chung Chang and Chih-Jen Lin. IJCNN 2001 challenge:
Generalization ability and text decoding. In Proceedings of International Joinl Conference on Neural Nelworks, volume 2, pp. 1031 1036. IEEE, 2001.
Christopher C Chapman, Didier P Monselesan, James S Risbey, Ming Feng, and Bernadette M Sloyan. A large-scale view of marine heatwaves revealed by archetype analysis. Nature Communications, 13(1):7843, 2022.
Yuansi Chen, Julien Mairal, and Zaid Harchaoui. Fast and robust archetypal analysis for representation learning. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1478 1485, 2014.
Adele Cutler and Leo Breiman. Archetypal analysis. Technomelrics, 36(4):338 347, 1994.
Adele Cutler and Emily Stone. Moving archetypes. Physica D: Nonkinear Phenomena, 107(1):1 16, 1997.
Anil Damle and Yuckai Sun. A geometric approach to archetypal analysis and nonnegative matrix factor- ization. Technometrics, (93):361-3700 2017.
Pierre De Handschutter, Nicolas Gillis, Arnaud Vandaele, and Xavier Siebert. Near-convex archetypal analysis. IEEE Signal Processing Letters, 27:81-85, 2019.
Michaäl Defferaad, Kirell Benzi, Pierre Vandergheynst, and Xavier Bresson. FMA: A dataset for music analysis. In 18th International Society for Music Information Retrieval Conference (ISMIR), 2017.
Maria R D'Esposito, Francesco Palumbo, and Giancarlo Ragozini. Interval archetypes: new tool for interval data analysis. Slatistical Analysis and Data Mining: The ASA Dala Science Journal, 5(4):322 335, 2012.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http:://arcive..isss... edu/
José H Dulá and Francisco J López. Competing nutput-sensitive frame algorithms. Computational Geometry, 45(4):186-195 197, 2012.
Irene Epifanio, M Victoria Ibáñez, and Amelia Simó. Archetypal analysis with missing data: see all samples by looking at few based on extreme profiles. The American Statistician, 74(2):169-183, 2020.
Flavia Esposito. A review on initialization methods for nonnegative matrix factorization: towards omics data experiments. Mathematics, 9(9):1006, 2021.
Manuel JA Eugster and Friedrich Leisch. Weighted and robust archetypal analysis. Computational Statistics & Dala Analysis, 55(3):1215-122) 2011.
Julia Mimbernat-Mayol, Albert Dominguez Mantes, Carlos D Bustamante, Daniel Mas Montserrat, and Alexander G Toannidis. Archetypal analysis for population genetics. PLoS Computational Biology, 18(8):
e1010301, 2022.
Donald Goldfarb and Shucheng Liu. An O(n primal interior point algorithm for convex quadratic pro- gramming. Mathematical Progrumming, 49(1):325 340, 1990.
Teofilo F Gonzalez. Clustering minimize the maxinum intercluster distance. Theoretical Computer
Science, 306, 1985.
13

---

Published in Transactions on Machine Learning Research (04/2024)
Ruijian Han, Braxton Osting, Dong Wang, and Yiming Xu. Probabilistic methods for approximate archetypal analysis. Information and Inference: A Journal of the IMA, 2022.
Abdelwaheb Hannachi and N Trendafilov. Archetypal analysis: Mining weather and climate extremes. Journal of Climate, 30(17):60276694 2017.
Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Courna- peau, Eric Wieser, Julian Taylor, Scbastian Berg, Nathanicl J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gerard-Marchann Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Ab- basi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 585(7825):
357-362, 2020.
Yuval Hart, Hila Sheftel, Jean Hausser, Pablo Szekely, Noa Bossel Ben-Moshe, Yael Korem, Avichai Tendler, Avraham E Mayo, and Uri Alon. Inferring biological tasks using pareto analysis of high-dimensional data. Nature Methods, 12(3):233- 225, 2015.
WK Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57 (1):97-109, 1970.
Jesper Love Hinrich, Sophia Elizabeth Bardenfleth, Rasmus Erbou Røge, Nathan William Churchill, Kristof- fer Hougaard Madsen, and Morten Mørup. Archetypal analysis for modeling multisubject fMRI data. IEEE Journal of Selected Topics Signal Processing, 10(7)11100-1171, 2016.
Dorit S Hochbaum and David В Shmoys. A best possible heuristic for the k-center problem. Mathematics of Operations Research, 10(2):100-184, 1985.
Catalin lonescu, Dragos Papava, Vlad Olaru, and Cristian Sminchisescu. Human3.6M: Large Scale Datasets and Predictive Methods for 3D Human Sensing in Natural Environments. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2014.
Leonard Kaufman and Peter Rousseeuw. Partitioning Around Medoids (Program PAM). Wiley And Sons, 1990.
Sebastian Mathias Keller, Maxim Samarin, Mario Wieser, and Volker Roth. Deep archetypal analysis. In German Conference on Pattern Recognition, pp. 171-185. Springer, 2019.
Sebastian Mathias Keller, Maxim Samarin, Fabricio Arend Torres, Mario Wieser, and Volker Roth. Learning extremal representations with deep archetypal analysis. International Journal of Computer Vision, 129 (4):805 820, 2021.
Yael Korem, Pablo Szekely, Yuval Hart, Hila Sheftel, Jean Hausser, Avi Mayo, Michael E Rothenberg, Tomer Kalisky, and Uri Alon. Geometry of the gene expression space of individual cells. PLoS Computational Biology, 117))ee0044224, 2015.
Laerke Gebser Krohne, Wang, Jesper Hinrich, Morten Mørup, Raymond CK Chan, and Kristoffer H. Madsen. Classification of social anhedonia using temporal and spatial network features from a social cognition fMRT task. Human Brain Mapping, 40(17):4965-4981, 2019.
Charles L Lawson and Richard J Hanson. Solving least squares problems, volume 15. SIAM, 1995.
Daniel D Lee and Scbastian Soung. Learning the parts of objects by non-negative matrix factorization.
Nature, (6755555) 791, 1999.
Stuart Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2):129-1377 1982.
Sebastian Mair and Ulf Brefeld. Coresets for archetypal analysis. In Advances in Neural Information Processing Systems, volume 32, 2019.
14

---

Published in Transactions on Machine Learning Research (04/2024)
Scbastian Mair, Ahcene Boubekki, and Ulf Brefeld. Frame-based data factorizations. In International Conference on Machine Learning, pp.
PMIRR, 2017.
Jieru Mci, Chunyu Wang, and Wenjun Zeng. Online dictionary learning for approximate archetypal analysis.
In Proceedings of the European Conference on Compuler Vision, pp. 486-501, 2018.
Morten Mørup and Lars Kai Hansen. Archetypal analysis for machine learning. In Proceedings of IEEE International Workshop on Machine Learning for Signal Processing, pp. 172-177. IEEE, 2010.
Morten Mørup and Lars Kai Hansen. Archetypal analysis for machine learning and data mining. Neurocom- puting, 80:54-63, 2012.
Susana Nascimento and Nuno Madaleno. Unsupervised initialization of archetypal analysis and proportional membership fuzzy clustering. In International Conference on Intelligent Data Engineering and Automated Learning, pp. 12 20. Springer, 2019.
Anders S Olsen, Rasmuus MT Høegh, Jesper L Hinrich, Kristoffer H Madsen, and Morten Mørup. Combining electro-and magneeeencceppallgaaphy data using directional archelypal analysis. Fronliers in Newroscience, 16, 2022.
Rafail Ostrovsky, Yuval Rabani, Leonard J Schulman, and Chaitanya Swamy. The effectiveness of lloyd-type methods for the k-means problem. Journal of the ACM, 59(6):1 22, 2013.
R. Kelley Pace and Ronald Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(3): 291 297, 1997.
Genevieve Patterson and James Hays. Sun attribute databasc: Discovering, annotating, and recognizing scene attributes. In IEEE Conference on Computer Vision and Pattern Recognition, pp. 2751-2758. TERE, 2012.
Sohan Seth and Manuel JA Eugster. Archetypal analysis for nominal observations. IEEE Transactions Pallern Analysis and Machine Intelligence, 38(5):849 861, 2015.
Sohan Seth and Manuel JA Eugster. Probabilistic archetypal analysis. Machine Learning, 102(1):85 113, 2016.
Jens Sjölund and Maria Bänkestad. Graph-based neural acceleration for nonnegative matrix factorization. arXiv preprint arXiv:2202. 00264, 2022.
Abdul Suleman. On illcooneiveed initialization in archetypal analysis. Advances in Data Analysis and Classification, 11(4):785 808, 2017.
Juliane Charlotte Thøgersen, Morten Mørup, Soren Damkiær, Søren Molin, and Lars Jelsbak. Archetypal analysis of diverse pseudomonas aernginosstranscriitomms reveals adaptation in cystic fibrosis airways. BMC Bioinformatics, 14(1):1 15, 2013.
Madcleine Udell, Corinne Horn, Reza. Zadch, and Stephen Boyd. Generalized low rank models. Foundations and Trends in Machine Learning, 9(1):1 118, 2016.
Andrew V Uzilov, Joshua M Keegan, and David H Mathews. Detection of non-coding RNAs the basis of predicted sccondary structure formation free energy change. BMC Bioinformatics, 7(1):1-30, 2006.
David van Dijk, Daniel B Burkhardt, Matthew Amodio, Alexander Tong, Guy Wolf, and Smita Krish- naswamy. Finding archetypal spaces using noural networks. In Proceedings of IEEE International Con- ference on Big Data, pp. 2634 2643. IEEE, 2019.
Guillermo Vinué, Irene Epifanio, and Sandra Alemany, Archetypoids: A new approach to define representa- tive archetypal data. Computational Statistics & Data Analysis, 87:102-115, 2015.
I-C Yeh. Modeling of strength of high-performance concrete using artificial neural networks. Cement and Concrete Research, 28(12):1797 1808, 1998.
Günter M Ziegler. Lectures on Polytopes. Graduate Texts in Mathematics. Springer New York, 2012.
15

---

Published in Transactions on Machine Learning Research (04/2024)
A Proofs
A.1 Proof of Lemma 3.1
Lemma 3.1. Let ||X - AZIl > 0, i.e., there are points gielding projection errors. Then, adding a point x E X \ Z to the set of archelypes Z according to AA++ (Algorilhm 7 is guaranleed to decrease the objective function.
Proof. Let P = (2)v(Z) be the polytope corresponding to the convex hull of the archetypes. There are only three cases. In the first case, the added point x lies within P. Then, P remains unchanged, and so do the projections of the points outside of P. Hence, the value of the objective function remains unchanged. However, AA++ (Algorithm 1) assigns a probability of zero such point x. Due to || AZl >> O, there has to be x E X' that has projection error and thus reduces the objective function when picked. In the second case, the added point x lies outside of P. and x is the only point projected on its face of P. Then, the objective function decreases exactly by the projection of x onto P. In the third case, the added point x lies outside of P and there are other points that lie on the same face of P as x. Then, adding x decreases the objective function by the projection of x itself and also the projection of the other points that lie on the same face. Overall, using AA++ (Algorithm 1) decreases the objective function with probability one.
A.2 Proof of Proposition 3.2
Proposition 3.2. Let / 21 1 and ||X AZAI > O. Per ileration, AA++ (Algorithm 7 gields a reduction
of at least as much as the reduction for the uniform énitialization in expectation. The reduction is equal only if every point in X \ Z has the same projection onto the current set of archetypes Z. Furthermore, using uniform initialization cannot reduce the objective function more than AA++ in expectation.
Within the proof we use Chebyshev's sum inequality which states that if p1 P 2 ... Pa and 2 2 Tri then Σn
E p.t; ( Ei ( 1 7T).
3.2 Idea: The cost that every point x, contributes to the objective function is d(x;, 2) = mingeoonv ||x; ql.. Hence, by adding x; to Z, the new objective function is reduced by at least d(x:, 2). If x; is the only point on its face of Z, the reduction only involves x, and hence the reduction is exactly d(x,, Z). If there are more points on the same face as X,, choosing x; reduces not only d(x,, Z) but also the projection of any other point on that face. Hence, d(x,, Z) serves as a lower bound on the reduction of the objective function when choosing x;. We now show that the reduction for AA++ is larger in expectation than for the uniform initialization.
Proof of Proposition
Define C' = Ei= d(xi, Z) which is equal to the normalization term within AA++. Without loss of generality, we can reorder the points such that d(x1, 2) 2 d(x2, 2) 2 ... 2 d(x,, 2). Let r; = d(x,, 2) be the cost (potential reduction) the point x; contributes to the objective function and p; = rr/C be the probability of choosing the point x; according to AA++ Thus, 71 2 r2 2 ... 2 Tn and p1 2 P2 2 ... 2 Pa. Furthermore, let u; = be the probability of choosing the point according to the uniform initialization. Due to Chebyshev's sum inequality, we obtain
1 i=1
piri
Pi
Pi
uet:
Ep[d(x;, 2)] E, [d( (x;, 2)),
i=1
by using E P. and by multiplying both sides by Hence, choosing x; according to AA++ (using the distribution p) results in a reduction which is at least large for a uniform initialization (using the distribution u). Equality is met, when 71 = f2 = ... Tu, hence
P2 = = Pn, and ultimately pt = for all i E [n). If any point has a larger distance than the others, we have rı > r2 due to the ordering and hence P1 > P2 in which the equality no longer holds and the reduction of the objective function is strictly larger for AA++ than for a uniform initialization expectation. Note that due to the construction of r, and Pi, we can never have the case of flipping the incquality. In other words, AA++ cannot reduce the objective function less than that of uniform initialization in expectation.
16

---

Published in Transactions on Machine Learning Research (04/2024)
B
Algorithms
Within the main body of the paper, omitted algorithmic pseudocodes whenever the procedure is not important or clear from the lext. However, we now provide those algorithmie descriptions for completeness. Algorithm 3 outlines the k-means++ initialization procedure, and Algorithm 4 shows the standard alternating optimization scheme of archetypal analysis according to Cutler & Breiman 1994).
Algorithm 3 +-means++ Initialization
1: Input: Set of n data points X, number of clusters k
2: Output: Initial set of & clusters centers 2
3: Sample index i uniformly at random from [n], i.e,, using p(i) n
4: Append X, to Z while [C] < do Sample i using p(i) x mingez ||x; - 2
Append x; to Z
end while
Algorithm 4 Archetypal Analysis (Cutler & Breiman, 1994)
1: Input: data matrix X RRdd number of archetypes k
2: Output: weight matrices A E Raxk and B E Rx* and archetypes Z E Rxd
3: Initialization of the archetypes E Rexd e.g., via AA++
4: while not converged do
a; arg min
IZ a; x,
Vie [n]
a11=1, a:20
6:
Z solve A ATAZ = A
2
bj arg min
| bj- 2
Vj [k]
b11=1, b,20
Z BX
end while
C Implementation Details
The code is implemented in Python using numpy Harris al. 2020 and is publicly available
For the optimization problem within AA++ and archetypal analysis, i.e., line 6 in Algorithm 1 and lines 5 and 7 in Algorithm 4 respectively, any QP solver can be used. In our implementation, we utilize the non- negative least-squares (NNLS) method from Lawson & Hanson (1995) and enforce the sunmation constraint by adding another equation the lincar system. For example, consider line in Algorithm 4 There, we add row of ones to ZT and one to X,, i.e., _ a 1-a; +1-00----- +2 .... 1- aa up one. Specifically, that the other lines of the linear system Z' a; x; do not dominate the summation to constraint. Note that using NNLS with a large constant M was already suggested by we use the Fortran implementation of NNLS from scipy i.e., @CCulelrr & Breiman (1994). Specifically, nnls, and slightly adapt it by increasing the hard coded maximum number of iterations since NNLS occasionally ran out of iterations.
ensure that the vector a sums scale the ones by M 1000 to have 1000 a1 1-000000ap 1:00-a2t.... 1000 -a 1000 such
D Results on Additional Data Sets
We further evaluate the initialization methods the following eight data sets. Airfoil (Brooks al 1989
has 1,303 data points represented in d = 5 dimensions. Concrete (Ych] 1998) has = 1, 030 instances in
httpss//scipy..rr
17

---

Published in Transactions on Machine Learning Research (04/2024)
Table 1: An overview of the 15 data sets used in this paper. The upper part is considered in the main body of the paper and the lower part is discussed in the appendix.
Data Set Name
Number of Data Points
Number of Dimensions
California Housing
20,640
Main Paper
Covertype
581,012
54
KDD-Protein
145,751
74
Pose
35,832
48
RNA
488,5
8
Song
515,345
90
FMA
106,574
518
Airfoil
1,503
Concrete
1,030
8
Appendix
Bankingl
4,971
7
Banking2
12,456
Banking3
19,939
Tjennl
49,990
22
MiniBooNE
130,064
50
SUN Attribute
14,340
102
d dimensions. The data sets Bonking1, Banking2, and Banking3 (Dulá López 2012] have n = 4, 971,
12, 456, and 19, 939 points in d 7, d = 8, and d = 11 dimensions, respectively. The Ijenn] data set has 49, 990 points in d 22 dimensions and was used in the LJCNN 2001 neural network competition [ We employ the same pre-processing as CCaagg & Lin (2001). MéniBooNE (Dua & Graff, 2017) consists of
= 130, 064 data points in d = 50 dimensions. Finally, the SUN Attribute (Patterson & Hays 2012) has
n = 14, 340 data points in d = 102 dimensions. A summary of all used data sets is provided in Table
Figure 7 depicts the performance for the additional eight data sets. Note that we omit AA++MC 1% for small data sets, i.e., if < 25, 000. Again, we can see that Uniform and the FurthestSum initializations often perform worse than the proposed approach and its approximation. AA+ almost always consistently best, except on some rare occasions, such as for k 25 on Ijennl. While the MCMC approximations of AA++ often perform very close to AA++ itself, we can some more significant performance gaps some of these additional data sets. Especially on BankProbleml, most Monte Carlo versions fail to approximate AA+ properly. However, note that all approximations are still better than the Uniform basclinc.
Figure summarizes aggregated statistics over 15 data sets, similar to Figure The difference that we now analyze the worst instead of the best-case setting. After initialization, the most frequently used initialization strategics for archetypal analysis, namely Uniform and FurthestSum, often yield the highest errors and, thus, the worst initializations. Overall, they sometimes recover from the sub-optimal initializations but still perform worse more often than the proposed AA++ strategy.
E Timing Results
The results depicted in Figures and 7 show the MSE as a function of iterations. However, as shown in Figure 0 the initialization times vary across all initialization strategies, and we should not assume that the time of an AA iteration is static, as discussed in Section G Thus, we also report on the MSE as a function of time in Figures 10 and 11 Surprisingly, although AA++ takes longer initialize, the overall optimization time for 30 AA iterations smaller for AA++ than for Uniform on Covertype, Concrete, and MiniBooNE data. A similar scenario can be seen for small & values on Song data. On some other data sets, using AA++ takes more overall, mainly due to the slower initialization. However, in turn, we get better performance. Furthermore, using an approximation of AA++ yields reasonable compromisc time-wisc while still being close to AA++ performance--iii...
hhttps://www.. ntu, wet.t.edduttw//cjiieee http:://wwccccaaaattttoooo/////a
18

---

Published in Transactions on Machine Learning Research (04/2024)
Usfor
Fist
FarthoASam
Morrod
M11
AA IMC 25%
MI IMC 10%
A4 IMC 5%
k-15
Airfoil 1-25
Airfoil 6-50
Airfoll A-75
1-100
10:
10
00-1
10.
101
10"
10
30
BonkProbiemi E-15
BankProslemi 2-25
BankProblem1 8-50
BankProblemi k-75
BoniProblemi A=100
10
10m1
10
00
10-1
101
10
10
10
10
BonkProblem2 A-15
BankProolem2 k-25
BankProblem2 8050
BankProblem2 k-75
BoniProblem2 1=000
10-1
10
10
1013
30
11
BankProblem3 15
BankProolem3 k--25
BankProblem3 500
BankProblem3 k5
BaniProblem3 A=100
10
101"
10
111"
111
11"
10
30:5
Concrete 15
Concrete 2525
Concrete 2=50
Concrete k5
Concrete A=
10
10
MSE
30
10
jonnl k-15
Jienl 2-25
Jinnl
Ijennl 75
[eennl 1=100
3x 10
3x 10"9
Sic uP1
M0E
4x 10
2x nr!
3xx 10
G 10
MiniBooNE h-15
MiniBooNE 2-25
MiniBooNE A-50
MiniBooNE k-75
Mi MiniBooNE 1-100
1.3
11F
10:
10
10
10
30-1
10-7
10"7
10
20
SUN Attribute E-15
SUN Attribute 25
SUN Attribute k-5
SUN Actribute 1575
SUN Actribute i-100
11
28 10
2.3 km
2x 10
1.5 x 10
3xx10
20x 109
23 6m
21 10
10x 10
24x
2 10
18x 11
22x
10 60
1.8 10
2x mf
13 1330000
18x
3.6 x 007
10c
10
15
20
15
20
20
20
10
30
Tteaaioos. of A
Rexaaior of AA
becations A
Iterodiess AA
Itorations A..
Figure 7: Results on Airfoil, Concrete, Bankingl, Banking2, Banking3, Ijennl, MiniBooNE, and SUN At- tribute using the CenterAndMarSeen pre- processcs as in the main body of the paper.
19

---

Published in Transactions on Machine Learning Research (04/2024)
Worst Initialization
Worst Overall
Median Initialization
Median Overall
A=
15
9
0
5
4
4
0
5
1
2
1
1
14
12
25
6
4
0
1
2
10
k=
0
3
4
4
1
1
0
k=75
11
0
4
0
2
6
2
0
4
1
0
3
3
0
100
9
0
5
0
0
3
5
2
1
4
1
4
1
0
Uniform
FurthestFiref
SurthestSum
AAcoresat
k-Means-t
AA+
Unifform
Furthest
AAcoreset
k-Means- +-+
AA++
Uniform
FurthestFing
FurthestSum
AAcoreset
k-Means-+
AA+-
FurthestFire
Uniform
FurchestSum
AAcoreset
A-Means-+
+- AA+++
Figure 8: Aggregated statistics over 15 data sets (seven data sets from above and eight data sets from the appendix) using CenterAndMuzSoule as pre-processing. Each table shows how often each initialization method yields the worst result for various choices & under different settings. Worst refers to the highest error of a single seed and median refers to the median over many seeds. We report on the performance after initialization and overall during the optimization.
Convex Hull of X
-2
-2
Data Set X
-2
4
2
2
-2
Figure 9: Comparison between the two pre-processing approaches. Left: Center data and divide by maximum value. Middle: Original data set. Right: Standardized data set.
F Pre- Pre-processing
In the main body of the paper, we pre-processed the data by first centering the data set and then dividing it by the maximum value. Another frequently used pre-processing scheme is standardization. That is, per dimension, we subtract the mean and divide it by the standard deviation. Since both are linear transfor- mations, they do not change the membership of the points being on the border the convex hull (Ziegler) 2012). However, the former scheme maintains the shape of the data set in terms of its convex hull, whereas the latter scheme changes it. This is depicted in Figure g where the original data is shown in the middle, the pre-processing of the main body of the paper is shown on the left-hand side and a data standardization is shown on the right-hand side.
G Results on Standardized Data Sets
We also conduct the same set of experiments on all data sets using standardization as pre-processing. In Figure 12] we can for the first set of data sets that FurthestSum usually performs worst, often by large margin. In contrast, the most consistent behavior has the proposed AA++ method, which is often the best. Besides, the proposed approximations of AA++ perform sufficiently close to AA++ itself. Figure 13 shows the same results but with respect to time instead of iterations. The performance of the second set of data sets is depicted in Figure 14 and the version showing time is provided Figure 15
Aggregated statistics on the results using standardization as pre-processing are depicted in Figure 16 As expected, AA++ wins on most of the data sets irrespective the applied setting. However, the numbers are slightly lower than in Figure 5 which summarizes the results using the Center AndMaScall pre-processing.
20

---

Published in Transactions on Machine Learning Research (04/2024)
Unioo
FasthoFFist
FattrsiSom
Murrss
AA1
AA MC 20%
AA MC 105
.... A.. MCC 5%
Califorii 2-15
Califoriia A-25
Calffrria 1-50
Califoriia A-75
California k-100
10
10
10"
10-1
10
10re
100
50
103
50
150
Computation Line soconds
Compocatón tine in scoonds
Corposattion Line sccurds
Comporatiinn line iccords
Computation time secords
Covertype 6-15
Covertype 2-25
Covertype i-50
Coveriype 1-75
Covertype Ao-100
101
10-1
117
10
10
10
10
20
25
50
100
100
200
500
Computatioo Eiee minutes
Computal on diee missee
Computation time IvEEUC
Computation time mino
Compucaton time miunes
FMA 2-15
FMA 25
FIM A-50
FMA A-
FMA K-100
101
A55E
10
Computatioon
Computation
Compucation
KDD Protein K-15
KDD Protein 2-25
KDD Protein R..50
KCD Protein A--75
KDD Protein A-100
11
100-2
10
MIE
10
10
10
21
40
Compatation Cre minetes
Competation Cinc in mEeuEe
Computation Sine minettes
Composatin lire minuee
Computation line 11 miluss
Pose K-15
Pose A--25
Pose k-50
Pose k-75
Pose 1-100
4
4.8 107
D. 100
3
3x
xx II1
3xx IIP"
10
5x 1
10
2xx 111
Six 100
25
07
100
100
200
100
300
300
470
500
Computtat om Tiee secoads
Computation on spo0ads
Computttion time secord
Compucation 1ine secorde
Computation time sccond
RNA K-15
RNA 2-25
RNA k-50
RNA A-
RNA 1o-100
111
10-4
10"
ASEE
10
1074
10
10
11
10
Computatioo
Compastaion
Song R-15
Song 2-25
Song A--5
Song h5-5
Song A-100
10
10
4x 10
8 10
110
201
1101
150
20
203 Compucation time
100
Comsutabion Emm minutes
Computat
Computation time
Computation time
Figure 10: Results on California Housing, Covertype, FMA, KDD-Protein, Pose, RNA, and Song, In contrast to F this figure shows the computation time of the initialization followed by 30 iterations of archetypal analysis on the a-axis.
21

---

Published in Transactions on Machine Learning Research (04/2024)
Unioom
Fast
FutinatSam
Aucroct
M-Marr
M4
AA: MC 201
AA: MC 10%
-... AA ICC 5%
k-15
Airfoil 2-25
Airfoil 1--50
Airfoi 6-75
Airfsi A-100
MSE
111
10
10
30
1011
10-1
10"
13-1
11777
10
10
0.0
015
1.0
2.5
2.
0
15
10
231
201
Computanicn on me seconds
Compusation time in secceds
Computabie time in secords
Computation time inconds
Computation time seceed
BankProbiemm 1=15
BandProbleml 2==2
BankProtiem] k=:50
BankProblem 7=
BankProšteml k=100
111
19:1
10
13
10
13
10
10
Comporation
Compozaation
Compotation
RenkPrabbemm 1-15
RankProblem2 J--5
RankProtiemO k-50
Pankeroblem22 2-75
PankProddeen 100
10
13
MSE
10
19
10-7
13+
10
20
100
t0
Computatoon on Line s0oodds
Corpozaioon line soccess
Computation Uree in socrdss
Computation line in secords
Compotation tine soconds
Barkeenneem3 2-15
BankProbilen3 2--25
BankPrabien3 k-50
Pankeroblem3 2-75
BankProtlen3 k-100
10
1012
10
13-2
10"
00
13
10
20
20
20
100
100
150
50
150
Compuration
Computatioo time
Computation seconds
Compuration time secords
A-25
Concrete K-50
Concrete R--55
Concrete K-100
117
10-2
MIE
10
10
In
0
10
10
20
Computse on time saconde
Compussion time secced
Computabor time
Computation time seconds
Computation time secodd
[ianel R-15
[icnn] 2-
[jann] A0
[icnni R-75
[jenn1 A-100
xx 20
xxx 10
10-1
4x 10
7 10
3
10-1-
100
20
100
100
100
200
200
100
Computat on socands
Compucation time
Compatation time secords
Computation time seconds
Compusation time accordh
MiniBooNE
MiniBooNE k=:55
MiniBooNE A=50
MiniBooNE As=75
MisiBooNE 1=:
13
10
10
30
10
10"5
10
10
13-1
Comsubationn minstee
Computation minute
Computation time in minstes
Compatation time minutee
Computation 1ime minuee
SUN Attribute k-15
SUN Attribute 2-25
SUN Attribute A-50
SUN Attrilute 4-75
SUN Attribute 100
2 10a
10
18x
28x 10
10
18% 1t.
1.8 109
22x 10
17 x 109
18% 11
x 10
11
82
L6 10
135 1 107
40 Computat on socands
30
101
200
130
300
D10
2001
320
Compusaiion time secoeds
Compatabinn secords
Computation time anconds
Computtaion time secodd
Figure 11: Results on Airfoil, Concrete, Banking1, Banking2, Banking3, Ijennl, MiniBooNE, and SUN Attribute using the CenterAndManSeule pre-processing as in the main body of the paper. In contrast to Figure 7 this figure shows the computation time of the initialization followed by 30 iterations of archetypal analysis on the axixs.
22

---

Published in Transactions on Machine Learning Research (04/2024)
FarthauFrnst
FarthoiSm
AA IMC
MC
Californis k=
Caiiforia 2=25
California A=
California k5=
Califoris k= k=100
-0"
117
107
N55
10
10"
10
Covertyps K-15
Covertype k-25
Covertype 1-50
Covertyoe k-
Covertype k-100
107
MSE
10P
-01
FMA 1=15
FMA 1-25
k-
FMA 1-100
107
110
ASEE
100
10
k-15
Proterr
KDD-Proter
107
10
6.
107
Sxx 101
3x 20/
2x 107
Pawe 1515
Pose
Pow
Pose J-
2
2x 101
3ix 100
10
NSE
6x 10f
101
8.x.0 U8
RNA 2515
RNA 2-25
RNA k-50
RNA A-
RNA K-100
'00
10"
10
10'
Song 6-15
Song 2-2
Song A-50
Song k-
Song 1-100
11P
ASE
102
107
10"
10'
Iiorations
locations
locutioro
Figure 12: Results on California Housing, Covertype, FMA, KDD-Protein, Pose, RNA, and Song using Standardization pre-processee pro-ssocsssing.
23

---

Published in Transactions on Machine Learning Research (04/2024)
Uniorm
FurthoatF ist
FurthacSom
Mcconeet
M-Maarai-
Mi-4
MA--MC
AA--MC
California 15
Califoenia 25
Calffrnia A=
California 7=
California A=100
10F
107
109
1888
10
10"
13
73
100
10C
100
Computation Binc accens
Computation lime accon
Compusation seconds
Compusation time seccess
Compotation time sccencs
Covertype Am115
Covertype 2-25
Covertppe Ž--50
Covertype A-75
Covertype 2=100
10a
MIE
11F-
11F
1000
Lin
Compotation
Computat
k-15
A--2
FMA A-500
FMA Jo-75
FMA k-1100
10f
107
M8E
107
10
20
0
50
50
70
50
300
150
100
203
300
100
200
300
400
Compuration time minues
Compusation time minuees
Comperation time in minurts
Compuraiin time in minotes
Computar on tire munes
KDD-Protein 15:15
KDD- Protein 2
KDD- Protein A=55
POD-Pratein 2::75
KDD-Protein Ro 100
112
6x 10
4x 12/
6x
10-
5x 12
x 10
8x 30
10
Compotation line minos
Computation time minus
Computation minuta
Compotatio Sanc ITTutO
Compotat Sirc
Posn 1-15
Poss 2-25
Poss k-50
Pose 6-75
Poss 6-100
2x 20f
MDE
2xx TIP
10
10
10
10
8x 10f
6x 107
100
200 Compozation
200 Compuction
Compusattion tine
lime
Compocation tine
RNA J-15
A-25
RNJ A-50
RNA A-75
A-100
31F
11f
M88
10
10
10
0.0
2.5
3.0
75
10.0
0
5
10
15
10
20
30
20
20
30
20
08
Compucatinn Time mvivutee
Compucation 1ime minevee
Compecation time in minures
Computaron vm time minruee
Computar Sire m 1170
Song 15:15
Song A5:25
Song 10
Song 75
Song A: 100
10
MEE
10P
11P
XC0
Compotation Line minatc
Compotation dinc minuss
Compotation in minata
Compotation tine minute
Compotation Einc IYDIO
Figure 13: Results on California Housing, Covertype, FMA, KDD-Protein, Pose, RNA, and Song using Stondardization pre-procss.nng. In contrast to Figure [2] this figure shows the computation time of the initialization followed by 30 iterations of archetypal analysis on the r-axis.
24

---

Published in Transactions on Machine Learning Research (04/2024)
FarthaeFrne
FustoosSa
Morroct
AA1-
Mi
IM-IMCAL
Airfai k==15
Airfll 2=25
Airfoil k=:50
Airfall Am=
Airfoil k= 100
1.2
10"
10-
10
10
11
10"2
100
RarkPPrhhemm1 1-15
RankPrcciemm k-25
Pankee 2-50
RankProbieem k-75
RanxProbleml 100
1.0
10"
1.0
107
10
1013
10"1
1073
10
101
107-2
107
10
RarkProhhem22 1-15
BankProisem? k-25
ParkProblem2 2-50
BankProliem2 k-75
PaniProbiem2 1-100
207
1.0
10r
10'2
10
Prrkees 1-15
BankPrdb k-25
ParkProbkm33 2-50
Praieen k-75
BaniProblem3 1-100
1.0
10F
1.0
20ml
10"
Concrete 1-15
Conceete A-25
Concrete 1-50
Concrete k-75
Concrete 1-100
2xx
10f
1.0
10
6xx
Scrrl k-15
Ijann] 2-
[jcnn1 k-50
[iamn] 1o
Ijenn] 2-100
111
107
1.2
15 11
x: 12
31
8.x10'' 10"
6x IIP
1x 22"
7 101
4x 10
1 107
10
MiniBooNE Am15
MinBBoNNE 2-25
MiniBooNE A-
Min BooNE k-
MiniBooNE k- 100
107
30
0.00
107
10-1
10
SUN Attribute 2-15
SUN Astribute 2-25
SUN Attribute k-50
SUN Astr A-
SUN Actribute 100
101
T 111
Dss 11P
1.0
10
Tx 10
6x 10
4 :00
4xx 10
6 101
ksL Iocatiors
118.
211
15
Rouuiios
Bicaliors
Figure 14: Results on Airfoil, Concrete, Bankingl, Banking2, Banking3, Ijennl, MiniBooNE, and SUN Attribute using Standordization pre-processing.
25

---

Published in Transactions on Machine Learning Research (04/2024)
Unior
FartheaFiist
FurthesSur
Mecroett
A-
AA1-MC
AA-+MC 10%
AA-MMC
Airfoil 1=15
Airfoi 25
Airooi 5=
Airfoil k=
Airfoil k=100
111
11-12
10"3
Compuration
Computation
Computaiion
Ranieechlemm] A-15
RankPrahiem] R-25
RankProlieml k-550
RandProbleml A--75
RankProbleml A-100
10
105
10"
10"1
10
10"5
10
1II
201
10
10
01
20
301
10
Compucaiion time seccace
Computatt on time in seconds
Comp ctmpitation time secords
Computation time in wcorch
Computanion tate In 180onds
Banidrroblem2 1::15
BankProblem22 2=
BankProblem2 k--50
BaniProblem2 15
BankProblem2 Ro 100
10
10
11
Compotation
RaniePrchlem3 1-15
RankProhiem3 k-25
RankProbiem3 k-50
RandProblem3 I=-55
PankProblem 1-100
101
1.0
10"
10
10
111
10
20
20
25
110
100
150
50
150
Comporation Sine SCC0ELS
Compatation on Uiee in socoads
Congulation Unic socord
Compusation line in acconds
Computai Lite in scc0mds
Concrete A-15
Concrete 2-25
Concrete k-50
Concrete A-75
Concrete 1-100
T0
11r
10
6xx. 10
10
10
10
10
15
Compucaion time
Compatation
Comp staiio time seconds
Computation thee seconds
Computation
Ijenn2
||can] A--2
Ijann1 A-50
Ijann] k-75
Ijenn1 k-100
107
10!
M3E
G 16f
4x 10
10f
2 Uf
50
100
170
Computation time vensce
Computaet on time secunds
Computiiicn time seconds
Computation time sccodd
Computaiion time tecoade
MiniBooNE 1-15
MiniBooNE A-25
MiniBooNE 8..50
MiniBoo k--75
MiniBooNE R-0100
10
11P
10
10
1.0
10"
10
10
10
10
010
25
5.1
T.b
1003
211
10
20
Computat on minue
Compatabon minutes
Compusttioon time minutes
Computaiion minutm
Computation time mirste
Attribute 15
Attribute k::25
SUN Astribute A=
Actribute 75
SUN Attribute i:
8
xxx. 001
7
38 9: 107
x 10
1.2
10"
10
6x 10
T> 100
10
7
bx 11F
6 10
25
100
1000 Computation time acconds
100
Compusation line sccuncc
Computation on socords
Computation Lime socond
Computatio in sccords
Figure 15: Results on Airfoil, Concrete, Bankingl, Banking2, Banking3, Ijennl, MiniBooNE, and SUN Attribute using Standardization as pre- pre-processing. In contrast to Figure 14 this figure shows the computation time of the initialization followed by 30 iterations of archetypal analysis on the z-axis.
26

---

Published in Transactions on Machine Learning Research (04/2024)
Best Initialization
Best Overall
Median Initialization
Median Overall
k=15
3
0
2
8
1
4
3
1
1
2
0
2
2
2
6
14
12
2=
3
0
2
3
0
3
2
7
110
5=
0
13
2
3
8
1
0
1
1
10
8
6
k=
75
0
1
0
1
1
12
2
0
1
2
10
1
0
1
0
13
2
0
0
2
1
10
100
0
2
11
3
0
3
8
0
0
13
2
0
0
0
0
13
2
FurthestFirst
Uniform
FurthestSum
AAcoreset
A-Means+++ +-+
AA++
Uniform
FurthestFirg
FurthestSum
AAcoreset
k-Means
+-+
AA++
Uniform
FurthestFirst
AAcoreset
k-Means+
+-
AA++
Uniform
AAcoreset
k-Means
+-+
AA++
Figure 16: Aggregated statistics over data sets (seven data sets from the main paper and eight data sets from the appendix) using Standardization as pre-processing. Each table shows how often each initialization method yields the best result for various choices of k under different settings. Best refers to the lowest single seed and median refers to the median over many seeds. We report the performance after initialization and overall during the optimization.
Worst Initialization
Worst Overall
Median Initialization
Median Overall
k=
0
9
3
3
0
1
1
8
1
5
4
0
0
14
12
-=
25-
0
11
4
2
0
0
8
0
5
3
0
0
100
50
9
4
0
7
6
8
6
k=
75
0
9
0
7
4
0
0
0
1
7
0
6
5
0
0
0
k= 100
8
5
3
0
0
0
1
7
0
7
4
4
0
2
Uniform FurthestSum
FurthestFirst
Accoreset
A-Means-++ +-+
AA+-+
Uniform
Accorset
k-Means++ +-+
AA+-
FurthestFire
Uniform
Furthest
AAcoreset
L-Means+ +-
AA++
Uniform
AAcoreset
B-Means++
AA++
Figure 17: Aggregated statistics over 15 data sets (seven data sets from the main paper and eight data sets from the appendix) using Standardization as pre-processing. Each table shows how often each initialization method yields the worst result for various choices of k under different settings. Worst refers to the highest error of a single seed and median refers to the median over many seeds. We report on the performance after initialization and overall during the optimization.
27


Reviews:
[{"review_id": "zPaipWXvY4", "reviewer": "Reviewer_fwRK", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. The paper theoretically shows that AA++ outperforms the uniform initialization.\n2. The empirical results show that the proposed initialization can help AA to converge to a better minima in many cases.\n\nWeaknesses:\n\n1. The theoretical contribution is a bit weak. Proposition 3.2 only says that AA++ is better than uniform sampling in expectation. But the paper does not describe how well the points sampled by AA++ (and AA++MC) are. Also, it's unknown whether the AA++ and AA++MC are better than other non-uniform initialization in theory. \n2. In figure 9 and 10, the paper shows that overall optimization time for 30 AA iterations is smaller for AA++ than for uniform initialization. However, the authors do not provide any explanation for this phnomenon.\n3. The dimension of the dataset is a bit small. I suggest the authors to perform same experiments on datasets with larger dimension (d>=500).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a novel initialization approach, called AA++, for archetypal analysis. The proposed method leverages the idea of kmeans++ and has better performance than uniform initialization. The authors also propose two strategies to approximate AA++ and significantly reduce the computational costs.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. For the theory, I would like to see more elaborate analysis for AA++ and AA++MC, as that of k-means++ (Arthur & Vassilvitskii, 2007) and k-MC^2 (Bachem et al. (2016).\n2. For the experiments, I would like to see both time comparison and performance comparison on datasets with larger dimension (d>=500).\n\nThe first suggestion is critical to my recommendation for acceptance and the second suggestion can strengthen the work in my view.", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "nKFbgUV8Ol", "reviewer": "Reviewer_tPNN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n- The paper is very clearly written. Although I am not familiar with archetypal analysis, I can still follow the paper clearly.\n- The paper includes very comprehensive experimental results. Experimental results include the comparison with existing baseline initialization methods as measured by MSE, the comparison on time and the aggregated statistics on the number of times the proposed method yields the best results. The paper considers the effect of the preprocessing strategy and the number of archetypal. The experimental results seem to be convincing. \n\nWeakness:\n- The theoretical results are not strong. The theoretical analysis only shows that the proposed method is no worse than the uniform initialization in expectation.\n- The proposed initialization method is time consuming. Fortunately, this problem is relaxed as the paper introduces two approximate methods.\n- The paper mainly considers the initialization method. While good initialization can improve the performance, a general belief is that the later archetypal analysis may be more important. In the experimental results, the paper only uses Algorithm 4 as the archetypal analysis method, which was proposed 30 years ago. Suppose we consider a more advanced archetypal analysis method, can we still get similar results showing the advantage of the proposed initialization method?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper studies the initialization strategy for archetypal analysis, which is a matrix factorization method with many applications in data science. The paper introduces the archetypal analysis++ initialization method. The basic idea is to build a distribution per iteration, which gives more weight to those examples far away from the convex hull of the existing archetype. Then, sample the archetypal from this distribution and append it. To further relieve the computational cost, the paper further proposes two methods: one avoids the computation of the distance to the convex hull, and the second uses MCMC to relax the dependency on the number of examples. The paper gives comprehensive experiments to show the behavior of the proposed initialization method.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- It would be helpful if the authors can add discussions why they consider the specific archetypal analysis method in Algorithm 4. The authors mentioned that each iteration in Algorithm 4 is expensive, and may be more expensive than the cost of the initialization, which justifies the importance of initialization. Suppose we consider a more efficient archetypal analysis method with less computational cost per iteration. Then is the initialization still important in this case?\n- Section 5: \"The same gape\" should be \"The same gap\"\n- Section 6: \"as a heuristic\" should be \"as a heuristic algorithm\"\n- Algorithm 4 requires to compute the inverse of $A^TA$. Is this matrix always invertible?", "broader_impact_concerns": "I have no broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "QChCKl3usn", "reviewer": "Reviewer_n3Fp", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The paper is extremely well written and is pleasant to read. I learned a lot from reading it. The figures are very helpful. \n- The proposed method is natural and sound\n- It indeed seems to work better than all the other considered baselines\n- The theoretical results look correct\n- The experiments are instructive and clearly show the benefits of the proposed methods.\n\nWeaknesses \n\nThe main weakness of that paper, in my view, is its theory. I find the theoretical part underwhelming. The authors propose an initialization scheme quite similar in spirit to k-means ++, yet they do not discuss similar theoretical results as the ones obtained in the original kmeans++ paper. The reader is left expecting a mention on the expected value of the loss function after initialization, as done in kmeans++. If the extension of the proof is hard and beyond the scope of the paper, then I think it should be mentioned in the paper, as well as the precise difficulties encountered by the authors.\n\nI have some other minor concerns:\n\n- In the experiments, some curves go up, doesn't the AA algorithm have a descent guarantee?\n- The algorithm used to solve the least-squares under simplex constraint could be explained in further detail, as it is a critical part of both AA and AA++ and is -- if I'm not mistaken -- the main bottleneck of the methods.\n\nMisc:\n- In fig2, recalling the MSE definition or at least pointing to it would help the reader\n- I think that fig 9 gives a more interesting picture than fig 4, since it clearly shows both the superiority of the proposed method and the advantage of the MC methods.\n- A parallel could be made between FurtherFirst and the proposed method in the following two ways. AA++ can be seen as a relaxation of furtherfirst, where instead of sampling from the argmax of a distribution the authors sample from the distribution itself. It can also be connected to a different loss function; indeed, in kmeans++, the authors mention that to optimize the loss $\\|X - ABX\\|^l_2$ then, one should sample from $ \\|X - ABX\\|^l_2$. If I'm not mistaken, FurtherFirst therefore corresponds to the $l\\to \\infty$ limit, being the optimal initialization for a cost function which is not the MSE.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a novel initialization scheme for AA (Archetypal Analysis), called AA++. \n\nGiven a dataset (x_1, \\dots, x_n), AA seeks $k$ anchor points $z_1, \\dots z_k$ that are each a convex combination of the $x_i$, so that each $x_i$ can be well approximated by a convex combination of the $z_j$. The loss function to minimize is $\\min_{A, B} \\|X - ABX\\|^2$ where $A$ and $B$ are entry-wise positive and row-wise stochastic.\n\nThis is a non-convex problem so initialization is of extreme importance. AA++ works in a similar fashion as k-means ++: the first point $z_1$ is taken uniformly among the dataset points, and then samples iteratively $z_{j}$ following a distribution proportionnal to $p(i)\\propto  \\min_{q\\in\\mathrm{Conv}(z_j)}\\|x_i - q\\|^2$.\n\nComputing each $p(i)$ is a quadratic program, hence the authors propose a cheaper version of the algorithm, by sampling only from a subset of points at each iteration.\n\nThe authors then compare their approach to several baselines on 14 datasets, and demonstrate the promises of their approach.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "In my view, either improving the theory of the paper or explaining clearly what is the limiting factor in extending the kmeans++ bounds on the loss would enhance the paper.", "broader_impact_concerns": "Na", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers praise the clarity and significance of the paper, and there is a clear unanimity to accept the paper. Congratulations on a fine piece of work!
Decision: Accept (poster)