extfieldsofstudy:
Geography
Computer Science

text:
Geometric methods for estimating representative sidewalk widths applied to Vienna’s streetscape surfaces database

Space, and in particular public space for movement and leisure, is a valuable and scarce resource, especially in today’s growing urban centres. The distribution and absolute amount of urban space—especially the provision of sufficient pedestrian areas, such as sidewalks—is considered crucial for shaping living and mobility options as well as transport choices. Ubiquitous urban data collection and today’s IT capabilities offer new possibilities for providing a relation-preserving overview and for keeping track of infrastructure changes. This paper presents three novel methods for estimating representative sidewalk widths and applies them to the official Viennese streetscape surface database. The first two methods use individual pedestrian area polygons and their geometrical representations of minimum circumscribing and maximum inscribing circles to derive a representative width of these individual surfaces. The third method utilizes aggregated pedestrian areas within the buffered street axis and results in a representative width for the corresponding road axis segment. Results are displayed as city-wide means in a 500 by 500 m grid and spatial autocorrelation based on Moran’s I is studied. We also compare the results between methods as well as to previous research, existing databases and guideline requirements on sidewalk widths. Finally, we discuss possible applications of these methods for monitoring and regression analysis and suggest future methodological improvements for increased accuracy.

Introduction
Today's IT capabilities and data collection on various urban management levels provide mobility relevant data on a comprehensive basis. In the past, random tape measurements provided anecdotal evidence of characteristic dimensions of transport infrastructure. Today's new geographic datasets enable systematic and comprehensive derivation of characteristic measurements. Besides detailed measurements for planning, the need of area-wide, quick, workable and relation-preserving overview methods may prove more important to some management needs than in-depth analysis. GIS is increasingly being used for deriving transport-related state variables, e.g. catchment areas of public transport stations (Iseki and Tingstrom 2014), the accessibility to public transport (O'Sullivan et al. 2000) or network performance (Mesbah et al. 2012). Walkability is another parameter, which has gained importance in recent quantitative analyses, for example analysing the connection between quality of urban built environment and mental health in youth (Duncan et al. 2013), active school transport (Wong et al. 2011), obesity (Duncan et al. 2012;Agampatian 2014) or the dichotomy of objective and perceived walking times (Dewulf et al. 2012). Findings suggest strong correlations between walkability and likelihood of taking a walking trip, vehicle miles travelled and obesity prevalence (Frank et al. 2007) and with moderateto-vigorous physical activity (Sallis et al. 2009). Due to the lack of more detailed data, a lot of studies stay on a meta-level when linking walkability with indicators such as landuse mix, residential density and intersection density (Frank et al. 2005) or pedestrian route directness. These can be calculated easily using GIS software and existing data. However, even simple indicators such as mere sidewalk presence or sidewalk widthused in a few studies (Ewing et al. 2004;Lin and Chang 2010)-need to be collected at least in part manually, which is time-consuming and costly (Schneider et al. 2005;Frackelton et al. 2013). Systemic data collection on sidewalk widths is needed. Here, our paper steps in with the proposition of automated computation of representative sidewalk widths with differently grained methods according to data availability.
Our paper is organized as follows: the next section presents three methods to compute representative sidewalk width and applies them to a dataset from Vienna, Austria. Section 3 presents distribution parameters of resulting representative widths and the mapped study of autocorrelation in terms of position in the city. The second to last section discusses the plausibility of results and draws conclusions on result significance elaborately. In the final section, results are reviewed as well as method application and improvement is debated.

The dataset
We use the SIS_F urban surface management database of the Vienna's road administration, which includes all public space surfaces within Vienna's municipal borders classified by type, pavement, maintenance responsibility, etc. This database contains 804,958 surface polygons and includes surface types such as road lanes, stairs, tramway rights of way, street islands and cycle tracks. Surface types that cover sidewalk areas are: sidewalk (GG), driveway entrance for house access (EE) and walkway entrance for house access (HH), which we use for our analysis. With few exceptions, these three types cover all walkable surfaces between facades and road lanes. The analysis does not include pedestrian precincts (FZ) since they tend to cover the total street width from facade to facade and separated walkways (FW) and since they denominate alignments detached from the regular streetscape, e.g. paths through parks. Table 1 describes the dataset (number of surfaces per type) as well as characteristic parameters (min, max, mean) of area and circumference.
In effect, minimum areas are not zero, but smaller than 100 cm 2 . As not only physical boundaries are shown in the SIS_F database, but also lines connecting house corners randomly delimit surfaces, such extremely small snippets exist. For example, the lower left corner of Fig. 2 shows such a triangular snippet-in some cases their size is less than the threshold given above.
A first analysis for Vienna's 23 administrative districts ( Fig. 1) suggests that increasing area size and distance from the city centre cause a sub-linearly growing portion of pedestrian surfaces, with an exponent of 0.77. This sub-linear development shows that peripheral districts with more recent road designs provide less pedestrian surfaces than central areas with often century old road layouts.

Three methods
We derived and applied three different methods for estimating representative sidewalk width w REP . Methods 1 and 2 calculate w REP per individual surface, while Method 3 aggregates on a street axis basis: • Method 1 uses the minimum circumscribed circle of a polygon.
• Method 2 uses the maximum inscribed circle of a polygon.
• Method 3 utilizes the total sidewalk area within a buffer around the street axis. The first two methods keep the distinction between the three surface types as given by the SIS_F database for reasons of precision and small-sized permanent changes in in situ widths. The different orientation of length versus width calls for the necessity to different methodical approaches. An inclusion of EE or HH surfaces to adjoining GG surfaces would reduce the detailed representation of constant width changes.
In Method 1, the representative sidewalk width is calculated as where A GG is the surface area of the type GG shape and d C is the diameter of the circumscribed circle. w REP,EE/HH is directly associated with d C as EE and HH type surfaces are mostly oriented across the logical sidewalk width. Since GG surfaces are mostly longitudinally oriented, w REP,GG is calculated from the area. To compute circumscribed circles, we use the maximum bounding geometry tool implemented in ArcGIS.
In Method 2, we use the maximum inscribed circle to calculate w REP as where A EE/HH is the surface area of the type EE or HH shape and d I is the diameter of the maximum inscribed circle. Contrary to Method 1, w REP,GG is directly associated with d I , while EE and HH surfaces are derived from the surfaces respective area. The centre of the largest circle that fits inside an arbitrary polygon is the point that is inside of the polygon and furthest from any point on the edges of the polygon. This largest circle touches the edges of the polygon on at least two points. The algorithm that we used is based on Garcia-Castellanos and Lombardo (2007) and was implemented in QGIS. The algorithm solves this problem iteratively to any arbitrary precision and ensures convergence towards a local maximum of the distance to the polygon edges, as follows: 1. Define a rectilinear search region R from (x MIN , y MIN ) to (x MAX , y MAX );  (2007) use N x = N y = 21; 3. Of the points inside the polygon, find the point that is furthest from any point on the edge; 4. From that point define a new R with smaller intervals and bounds, where Garcia-Castellanos and Lombardo (2007) reduce R by a factor of ffiffi ffi 2 p ; and repeat from step two to get to any arbitrary precision. Figure 2 illustrates the application of Methods 1 and 2 to type GG, HH and EE surfaces.
Rather than looking at each sidewalk area on its own, the buffered street axis approach of Method 3 evaluates the total combined sidewalk, walkway and driveway area for each street. The representative width is computed as Herein, A is the total sidewalk area (all surface types combined) within the buffer, l is the length of the street axis and b is the buffer size. For the evaluation in this paper, we use a buffer of 20 m around each street axis feature. As shown in Fig. 3, the buffers extend beyond the ends of the street axis to include the intersections with all of their sidewalk areas. This can lead to an overestimation of w REP , particularly for short streets. While 20 m are sufficient to enclose the sidewalks of most streets in Vienna, some broad boulevards-such as the Ringstraße-would require bigger buffers. Because of the absence of data about actual total street width, the analysis therefore has to carefully weigh between missing some sidewalks due to smaller buffer sizes or bigger distortions of w REP for short streets due to bigger buffer sizes.

Circumference over area
For the large scale sidewalk analysis, we adopt a scale-free approach and plot circumference over area of every surface type . We call them C-A diagrams in short. Regular shapes with the shape's circumference U being a function of area A, according to Eq. (4) facilitate the interpretation of the results and are therefore added to the C-A diagrams.
Table 2 provides the parameters for these regular shapes, representing a circle, a square, a rectangle with a side ratio of 1:2 and a rectangle with 1:10. The area's exponent b remains constant at 0.5, the square root obviously. Figure 4 a-c shows that sidewalk surfaces (GG) are larger than driveway entrances (EE) and walkway entrances (HH) with a tendency to grow narrower with increasing area. This is shown in Fig. 4a by the grey cloud and contour lines moving away from the circle line (red) over the rectangle 2a line (orange) to the rectangle 10a line (purple) with surface sizes around 10 2 m 2 and larger. Figure 4a shows for sidewalk surfaces (GG) that density isolines clearly group along the orange line with the peaks in an area range of 10 0.5 -10 1.5 m 2 , this is from 3 to ca. 32 m 2 . Figure 4b illustrates that driveway entrances (EE) are bigger in size than walkway entrances (HH, Fig. 4c), as one would expect from door widths. Density peaks in a range around 10 1 m 2 for walkway entrances versus a range between 10 0 and 10 1 m 2 for driveway entrances, which peak around 10 0.5 m 2 of area and a 10 1 m circumference. The density spikes clearly group along the rectangle 2a line with no grey circle actually touching the red circle line. Walkway entrance (HH) density groups between the circle (red) and the rectangle 2a line (orange) as well and it can be observed that no grey dot reaches as far as the circle line. Both surface types, EE and HH, are distinctly smaller than GG surfaces and appear to be packed around the rectangle 2a line, whereas higher densities (contour lines) of the GG surfaces reach and cross the rectangle 10a line. In comparison with GG surfaces, EE and HH surfaces show a more compact form, as their density peaks appear between the lines for circles and rectangles with a side proportion of one-two times a. Although, from this point of view, it appears likely to expect that the estimation method for types EE and HH may lead to greater deviations of w REP from the actual sidewalk width, we apply this method nevertheless since the C-A plots indicate that the three applied methods are legitimate and valid. On the one hand side, for GG surfaces, Method 1 is expected to underestimate w REP , because the rectangular-style polygon's area is divided by its diagonal, which is bigger than the longitudinal side. Method 2, on the other hand, is expected to overestimate widths because it relies on the maximum circle. For types EE and HH, the circumscribed method, due to seeking for the polygon's longest diagonal, is expected to overestimate w REP . Method 3 is expected to behave indifferently, as the buffer includes sidewalk surfaces from the intersecting roads and the buffer width is added to street length as denominator in Eq. (3).

Results
The mean resulting values and the 15th and 85th percentile of resulting w REP are shown in Table 3. Percentile values were chosen over extremal values (min, max) as the SIS_F database includes surfaces close to zero in area and circumference-see Table 1.
Method 3 aggregates from n = 519,034 surfaces to n = 26,639 street links. Table 4 summarizes the cumulative results by giving the spans for w REP,15% , w REP,50% and w REP,85% values distinguished by administrative districts. The 85th percentile was chosen, as this is a standard percentile in transport measurement and design. For example v 85 gives the speed that is underrun by 85% of drivers/riders, a parameter often used for road design (FSV 2011(FSV , 2014(FSV , 2015a. The widths exceeded by 85% of measurements are given by w REP,15% . In addition, the array of districts that yield the envelope curve of maximum and minimum w REP cumulative distribution lines is given in the last two columns. Vienna's district numberings runs from the inside out, from 1 to 23:1 is the central business district, 3, 5 and 6 are within the second ring road, while the city's fringe comprises of districts 13 and 14 in the West, 22 in the East and 23 in the South.
As the aggregation of w REP per city district proves to deliver limited further insight, we calculate means of w REP per cell of a 500 by 500 m grid over the city's extents. Figures 5 and 6 illustrate the spatial distribution of mean w REP for Methods 2 and 3, respectively. Method 1 is not depicted since the resulting visual is very similar to Method 2 as shown in Fig. 5.
Visual analysis of the distributions suggests a non-random distribution. We therefore check for spatial autocorrelation using Moran's I analysis (Li et al. 2007). The Moran's I value for method 3 is 0.67, while method 1 and 2 result in a Moran's I of 0.31 and 0.34, respectively. A Moran's I value of -1 indicates perfect dispersion, while ?1 indicates perfect correlation, with values close to -1/(n -1) indicating a random spatial pattern. 1 Both z-scores and p values of all three methods indicate spatial autocorrelation that is significant at the 5% level. For Method 3 (buffered street), a highly significant cluster of wide sidewalks emerges in the inner city, while this is not observable for Methods 1 and 2. Figure 7 shows hot and cold spots detected using a local Moran's I analysis of Method 3 results. Dark red cells indicate hotspots where sidewalks are wide in the respective cell as well as in its neighbouring cells. These cluster in the center of the city dark blue cells at the city Type GG surfaces tend to have a different logical orientation as elements of the sidewalk than EE and HH surfaces: GGs are oriented along the sidewalk, while EEs and HHs are usually oriented across. Therefore, we plot cumulative distributions of GG versus EE ? HH samples separately for Methods 1 and 2-circumscribed and inscribed (Fig. 8). The difference of Vienna's total sample (thick lines) at 50% shares is 1.50 m for circumscribed and 0.75 m for inscribed. Table 5 shows the width ranges at specific percentages for comparison between datasets and Methods 1 and 2. Method 2 appears to close in from Method 1, as ''black'' moves up the scale, while ''blue'' moves down. Method 3 could not be analysed in such a way by definition, since it does not distinguish between surface types.

Discussion
Results are as granular as the input data and method characteristics determine. In our results, we have aggregated all representative sidewalk widths, either by districts or in a grid. These are two options among many, and for other use cases it may be more appropriate to aggregate on more detailed levels, such as on a building block or street level or even individual surfaces.

Context
While the results of our methods agree on some aspects, considerable differences between different method results can be observed as well, which we have presented in graphical and numerical form. Because the cumulative distributions for Vienna (Fig. 8) take a sigmoidal shape, which is often an indicator for a normal distribution, we tested it. Both, the visual evaluation of histograms as well as the Kolmogorov-Smirnov normality test, which is considered suitable for datasets bigger than 2000 specimen, reveals that all three w REP data sets are not normally distributed. The cumulative distributions of w REP (Fig. 8) reveal a long tail reaching to 10 m of width and beyond. These surfaces do not necessarily represent classic sidewalk settings (as in a cross section of wall-sidewalk-lanes-sidewalk-wall) but represent squares and plazas that are not classified as type FZ (pedestrian precinct). The distributions' long tails are one reason why the normality tests do not reveal a Gaussian distribution of data. There is a significant difference in results when using Method 3 in comparison with Methods 1 and 2, where the circumscribed method results in the least spreading/diverging curves. While individual district results from Method 3 vary considerably from Methods 1 and 2, spanning a wider range of values, the cumulative curve for Vienna's total datasets does not differ much. District distributions show a relative stability explicated as envelope curves to the set of curves. While district 23, a peripheral district with a lot of industrial areas, consistently appears on the list of minimum sidewalk widths, other peripheral districts (such as 13, 14 and 22) are part of the smallest widths curve sets as well. Maximum widths show more variety. Nevertheless, districts 1 and 20 bound the curve sets from the upper end in two out of three methods, indicating a intermethodical stability. Figure 9a clearly shows that, for type GG, the inscribed circle method overestimates w REP in comparison with the circumscribed method. For types EE and HH (Fig. 9b), an opposite relation is observable, and the inscribed method underestimates in comparison with circumscribed. Both observations are in line with expectations, see scheme in Fig. 2. The linear correlation coefficients for EE and HH are smaller than 1 (0.83 and 0.84), while GG surfaces result in a larger than 1 correlation (1.17). All three correlations display very high goodness of fit values (R 2 ) with 0.844 (HH), 0.896 (EE) and 0.904 (GG). In contrast to GG surfaces, EE and HH surfaces display a rather strict superior boundary at the relation of w REP;I =w REP;C ¼ 1 with just a tiny proportion of outliers beyond this virtual border. The relation of circumscribed to inscribed for GG surfaces reveals a larger than 1 coefficient, but does not show any boundary, as the dots appear to scatter more symmetrically around the first median line.
In comparison with single surface consideration in Methods 1 and 2, the buffered street axis method (3) provides an aggregate approach. Surfaces from all three types within the buffer are consolidated and divided by the street buffer's characteristic length l ? 2b. For sake of comparison with Methods 1 and 2, we assume most streets to have a left and a right sidewalk; therefore, we multiply with 0.5. The use of l ? 2b instead of simply l as denominator tends to produce better results since otherwise w REP would be overestimated systematically. This means that circumscribed and inscribed methods favour streets with one wide sidewalk over those with two narrower sidewalks. An improvement of the buffered street axis method would ask for the variation of buffer size within a sensitivity analysis of results.
When street link datasets include lane widths and lane numbers, buffer size could be adapted accordingly for every individual link. This would improve the accuracy of results considerably, instead of running a constant buffer size analysis over every type and dimension of roads. To our knowledge, only two approaches related to ours exist in the literature. One is Freeman and Shapira's (1975) development of the minimum area rectangle (MAR) method for calculating the length and width of a rectangle circumscribing an entity. It is likely that this rectangle's length reveals similar results to our circumscribing circles-Method 1. In the other approach, Zhu and Lee (2008) investigated Moran's I for various transport-related variables. While some variables showed a significant impact, sidewalk width did not show significant clustering.

Comparison and validation
To further evaluate our results, we compare our methods to an independent reference dataset called graph integration platform (GIP). GIP is the public administration's official street network dataset of Austria. 2 Each street network link is described in detail, including cross-sectional information with street and sidewalk average and minimum width. The average sidewalk width along a street link l, which is used for the following comparison, is computed as the sum of sidewalk widths w k over all sidewalks along l-Eq. (5) where p denotes the percentage to and from and describes the exact location of the modelled sidewalk along the link. Figure 10 presents a scatterplot of w REP according to the buffered street axis method and the width according to GIP. Each data point represents one network link. The observed R 2 is 0.504. The comparison shows that GIP widths tend to exceed results of the buffered street axis method. The biggest errors are due to the following effects: (1) streets where the 20 m buffer only contains parts of the relevant sidewalk areas; (2) differences between GIP and the SIS_F urban surface management database, where GIP tends to be more generous while the SIS_F database tends to classify more areas as non-sidewalks (for example, green areas surrounded by sidewalks). So far, none of our methods consider the qualitative distinction between a wide one-sided sidewalk and narrow two-sided sidewalks. Even though the resulting numbers may be the same, e.g. in Method 3, the qualitative perception by pedestrians may differ considerably.
To get a feeling for the quality of our applied methods, we performed a hands-on plausibility check. For each surface in a selected neighbourhood (three blocks), we manually measured the sidewalk widths in the SIS_F data set. For nearly rectangular surfaces, this was the smallest distance perpendicular to the kerbstones. For irregular surfaces, the distance most plausibly perceived as representative sidewalk width was measured. The plausibility check delivered the following results (see Fig. 11): • Finding 1: w MEASUREDw REP differs most for GG; HH produces the best fit, EE slightly less. • Finding 2: EE surfaces (driveway entrances) are often wider than the sidewalk width. • Finding 3: The inscribed method is more exact than the circumscribed method.
Finding 1 is due to the fact that GG surfaced are often irregular, while HH and EE surfaces are mostly perfectly rectangular.
Finding 2 results in an overestimation of w REP as the surface orientation appears contrary to the method's preconditions-in partial contrast to the conclusion drawn in Sect. 2.3. We therefore calculated a corrected width ''EE*'' with the method used for GG rather than HH which led to a significant increase in fit. Finding 3 is a consequence of the used calculation method where especially in rectangular surfaces the inscribed method renders perfect results for w REP , while the circumscribed method overestimates distances due to the Pythagorean theorem. Given a sufficient number of HH and EE surfaces, this (adapted) method produces a good overview of w REP values.
One further approach for improvement of the presented methodology is to systematically calculate the logic transversal and longitudinal orientation of sidewalk surfaces from neighbourhood relationships between surface types irrespective of actual surface orientation. This is expected to be an advancement benefitting from treating groups of surfaces instead of every surface by itself. Under ideal circumstances, a GG type surface is expected to have its longer extension oriented with the longitude of the sidewalk, while for surfaces denoting the access to buildings over the sidewalk, such a clear distinction between width and length in terms of pedestrian's logic is not possible from geometry alone.
Sidewalk surfaces, especially in intersection setups, are segmented by lines that result from virtual alignment lines between house corners. These fragmentations reduce surface size artificially and are likely to lead-especially for Method 2-to an underestimation of w REP . For example, when a square is divided in two isosceles triangles, the relation of areas is 2, but the relation of inscribed diameters changes to 1.707. Circumscribed diameters remain the same.

Requirements
Measurements of the actual space that pedestrians demand for (Schopf 1985), put our results into the following perspective: two-people encounters ask for a width of 2.21 m at 50% share and 2.56 m at 85% share. While with the circumscribed method the 50% value is met only by the better half of districts, the 85% value is easily overachieved. A similar situation is present with the inscribed method. The buffered street axis method results in show that the total dataset meets the 50% as well as the 85% requirements. The Austrian transport engineering guidelines on sidewalk design (FSV 2015a, b) give 1.50 m as the minimum width for enabling the meeting of two opposite walking pedestrians and ask for a basic sidewalk width of 2.0 m. All three methods result in a w REP,50% that meets this requirement, indicating that only selected districts assessed using the inscribed method (w REP,15% up to 1.51 m), do not meet the quality standard for sidewalk infrastructure (Table 3).

Conclusion
In this paper, we have used three different methods for calculating a representative sidewalk width from a detailed city-wide urban surfaces GIS database.
Our paper contributes insofar as it introduces three methods which can be used to quickly estimate representative sidewalk widths from existing large scale urban surface GIS data. Moreover, the methods' quality, significance and representativeness are tested in a series of critical examinations. In one case, surface types EE, the comparison with manual measurements suggest to change the calculation algorithm accordingly when applying Methods 1 and 2 as these surfaces appear to be more oriented like GG surfaces than HH surfaces.
Results are a quantitative basis for subsequent in-depth analysis. Not only can these results be precisely compared to guideline requirements, they can also evaluate in terms of spatial autocorrelation. The question ''Do clusters of extra wide or extra narrow sidewalks exist?'' can be answered for areal entities. We performed this examination based on a 500 by 500 m grid to show that Method 3 leads to such a clustering of likewise results.
The task, which cannot be satisfied with our methods, but is one of very high importance to real pedestrian quality, is the consideration of net sidewalk widths. Gross width encompasses the sidewalk surface from wall to kerb, while net width is determined by obstacles and hindrances such as traffic sign posts, bargain baskets, light poles and advertisements cluttering sidewalks. As Fig. 2 shows in its left lower corner, the database already includes some of these fixed obstacles: tree bases and traffic light control boxes. Other fixed obstructions are not included in the present data base. To our knowledge, only lighting post GIS data exists in Vienna to a level precise enough to be included in our kind of calculations-a task open for future endeavours. One pure geometrical approach would be to cookie-cut these obstructions from its underlying sidewalk surfaces. Yet, most of these obstructions such as lamp poles and traffic signs engage only a very small horizontal space, hardly impacting on surface size. A lamp pole's obstruction effect lies in its local space constriction-a singularity in comparison with sidewalk surfaces. Although the road engineering rules pose a car-lane envelope driven demand for putting traffic signs up to 60 cm into the sidewalk surface, systematic surveys of how much these obstacles reduce gross widths have not been undertaken so far to our knowledge. But, we consider research on this issue as dearly needed. In a first approach, undergrad students at our Institute of Transportation are now conducting preliminary enquiries addressing net and gross sidewalk width discrepancy.
In general, GG surfaces due to their variety in shapes and sizes promise lower accuracy than EE and HH surfaces. But how representative are EE and HH surfaces-in their appearance discrete information-in contrast to the larger GG surfaces, representing most of the sidewalk surfaces? Therefore, the question of weighting of method's results in aggregate analysis needs to be targeted for method improvement. As of now, every derived w REP value weighs in equally. However, a weighting method for calculating w REP for road sections or building blocks has to be defined. Shall GG based results are forfeited for the more precise but likely less representative EE-and HH-based values? It appears conceivable to the authors that when using GG, EE and HH-based results, the individual surfaces size ought to be included as a weight in calculating aggregated representative values, e.g. for a block, a street section or an administrative unit. The question, when is it sufficient to just arithmetically average EE and HH values to achieve an acceptable degree of representation, poses a research task in need of being addressed. We see adding improved methods for dealing with all walkable surfaces-such as pedestrian precincts and separated walkways-as one additional horizon of improvement.
For the methods presented, we envision the following application examples: 1. Automatic enrichment of street network databases with sidewalk width information (for example for the official Austrian reference network GIP as an alternative to extensive manual data collection). 2. City-wide systemic monitoring of sidewalk widths-as absolute values or as a proportion of the total streetscape. 3. Input data for surveys, e.g. Frank et al. (2007).
Thus, we consider our three methods as quick tools for enabling a better understanding of urban systems at the crossroads of transport and surface allocation-thus a step towards leaving the yardstick behind.


extfieldsofstudy:
Physics
Mathematics
Computer Science
Biology
Economics
Medicine

text:
*K-means and cluster models for cancer signatures

We present *K-means clustering algorithm and source code by expanding statistical clustering methods applied in https://ssrn.com/abstract=2802753 to quantitative finance. *K-means is statistically deterministic without specifying initial centers, etc. We apply *K-means to extracting cancer signatures from genome data without using nonnegative matrix factorization (NMF). *K-means’ computational cost is a fraction of NMF’s. Using 1389 published samples for 14 cancer types, we find that 3 cancers (liver cancer, lung cancer and renal cell carcinoma) stand out and do not have cluster-like structures. Two clusters have especially high within-cluster correlations with 11 other cancers indicating common underlying structures. Our approach opens a novel avenue for studying such structures. *K-means is universal and can be applied in other fields. We discuss some potential applications in quantitative finance.

Introduction and summary
Every time we can learn something new about cancer, the motivation goes without saying. Cancer is different. Unlike other diseases, it is not caused by "mechanical" breakdowns, biochemical imbalances, etc. Instead, cancer occurs at the DNA level via somatic alterations in the genome structure. A common type of somatic mutations found in cancer is due to single nucleotide variations (SNVs) or alterations to single bases in the genome, which accumulate through the lifespan of the cancer via imperfect DNA replication during cell division or spontaneous cytosine deamination [1,2], or due to exposures to chemical insults or ultraviolet radiation [3,4], etc. These mutational processes leave a footprint in the cancer genome characterized by distinctive alteration patterns or mutational signatures.
If we can identify all underlying signatures, this could greatly facilitate progress in understanding the origins of cancer and its development. Therapeutically, if there are common underlying structures across different cancer types, then a therapeutic for one cancer type might be applicable to other cancers, which would be a great news. 2 However, it all boils down to the question of usefulness, i.e., is there a small enough number of cancer signatures underlying all (100+) known cancer types, or is this number too large to be meaningful or useful? Indeed, there are only 96 SNVs, 3 so we cannot have more than 96 signatures. 4 Even if the number of true underlying signatures is, say, of order 50, it is unclear whether they would be useful, especially within practical applications. On the other hand, if there are only a dozen or so underlying signatures, then we could hope for an order of magnitude simplification.
To identify mutational signatures, one analyzes SNV patterns in a cohort of DNA sequenced whole cancer genomes. The data is organized into a matrix G is , where the rows correspond to the N = 96 mutation categories, the columns correspond to d samples, and each element is a nonnegative occurrence count of a given mutation category in a given sample. Currently, the commonly accepted method for extracting cancer signatures from G is [5] is via nonnegative matrix factorization (NMF) [6,7]. Under NMF the matrix G is approximated via G ≈ W H, where W iA is an N × K matrix, H As is a K × d matrix, and both W and H are nonnegative. The appeal of NMF is its biologic interpretation whereby the K columns of the matrix W are interpreted as the weights with which the K cancer signatures contribute into the N = 96 mutation categories, and the columns of the matrix H are interpreted as the exposures to the K signatures in each sample. The price to pay for this is that NMF, which is an iterative procedure, is computationally costly and depending on the number of samples d it can take days or even weeks to run it. Furthermore, it does not automatically fix the number of signatures K, which must be either guessed or obtained via trial and error, thereby further adding to the computational cost. 5 Some of the aforesaid issues were recently addressed in [8], to wit: (i) by aggregating samples by cancer types, we can greatly improve stability and reduce the number of signatures; 6 (ii) by identifying and factoring out the somatic mutational noise, or the "overall" mode (this is the "de-noising" procedure of [8]), we can further greatly improve stability and, as a bonus, reduce computational cost; and (iii) the number of signatures can be fixed borrowing the methods from statistical risk models [9] in quantitative finance, by computing the effective rank (or eRank) [10] for the correlation matrix Ψ ij calculated across cancer types or samples (see below). All this yields substantial improvements [8].
In this paper we push this program to yet another level. The basic idea here is quite simple (but, as it turns out, nontrivial to implementsee below). We wish to apply clustering techniques to the problem of extracting cancer signatures. In fact, we argue in Section 2 that NMF is, to a degree, "clustering in disguise". This is for two main reasons. The prosaic reason is that NMF, being a nondeterministic algorithm, requires averaging over many local optima it produces. However, each run generally produces a weights matrix W iA with columns (i.e., signatures) not aligned with those in other runs. Aligning or matching the signatures across different runs (before averaging over them) is typically achieved via nondeterministic clustering such as k-means. So, not only is clustering utilized at some layer, the result, even after averaging, generally is both noisy 7 and nondeterministic! I.e., if this computationally costly procedure (which includes averaging) is run again and again on the same data, generally it will yield different looking cancer signatures every time!
The second, not-so-prosaic reason is that, while NMF generically does not produce exactly null weights, it does produce low weights, such that they are within error bars. For all practical purposes we might as well set such weights to zero. NMF requires nonnegative weights. However, we could as reasonably require that the weights should be, say, outside error bars (e.g., above one standard deviationthis would render the algorithm highly recursive and potentially unstable or computationally too costly) or above some minimum threshold (which would still further complicated as-is complicated NMF), or else the noncompliant weights are set to zero. As we increase this minimum threshold, the matrix W iA will start to have more and more zeros. It may not exactly have a binary cluster-like structure, but it may at least have some substructures that are cluster-like. It then begs the question: are there cluster-like (sub)structures present in W iA or, generally, in cancer signatures?
To answer this question, we can apply clustering methods directly to the matrix G is , or, more, precisely, to its de-noised version ′ G is (see below) [8]. The naïve, brute-force approach where one would simply cluster G is or ′ G is does not work for a variety of reasons, some being more nontrivial or subtle than others. Thus, e.g., as discussed in [8], the counts G is have skewed, long-tailed distributions and one should work with log-counts, or, more precisely, their de-noised versions. This applies to clustering as well. Further, following a discussion in [11] in the context of quantitative trading, it would be suboptimal to cluster denoised log-counts. Instead, it pays to cluster their normalized variants (see Section 2 hereof). However, taking care of such subtleties does not alleviate one big problem: nondeterminism! 8 If we run a vanilla nondeterministic algorithm such as k-means on the data however massaged with whatever bells and whistles, we will get random-looking disparate results every time we run k-means with no stability in sight. We need to address nondeterminism! Our solution to the problem is what we term *K-means. The idea behind *K-means, which essentially achieves determinism statistically, is simple. Suppose we have an N × d matrix X is , i.e., we have N dvectors X i . If we run k-means with the input number of clusters K but initially unspecified centers, every run will generally produce a new local optimum. *K-means reduces and in fact essentially eliminates this indeterminism via two levels. At level 1 it takes clusterings obtained via M independent runs or samplings. Each sampling produces a binary N × K matrix Ω iA , whose element equals 1 if X i belongs to the cluster labeled by A, and 0 otherwise. The aggregation algorithm and the source code therefor are given in [11]. This aggregationfor the same reasons as in NMF (see above)involves aligning clusters across the M runs, which is achieved via k-means, and so the result is nondeterministic. However, by aggregating a large number M of samplings, the degree of nondeterminism is greatly reduced. The "catch" is that sometimes this aggregation yields a clustering with K′ < K clusters, but this does not pose an issue. Thus, at level 2, we take a large number P of such aggregations (each based on M samplings). The occurrence counts of aggregated clusterings are not uniform but typically have a (sharply) peaked distribution around a few (or manageable) number of aggregated clusterings. So this way we can pinpoint the "ultimate" clustering, which is simply the aggregated clustering with the highest occurrence count. This is the gist of *K-means and it works well for genome data.
So, we apply *K-mean to the same genome data as in [8] consisting of 1389 (published) samples across 14 cancer types (see below). Our target number of clusters is 7, which was obtained in [8] using the eRank based algorithm (see above). We aggregated 1000 samplings into clusterings, and we constructed 150,000 such aggregated clusterings (i.e., we ran 150 million k-means instances). We indeed found the "ultimate" clustering with 7 clusters. Once the clustering is fixed, it turns out that within-cluster weights can be computed via linear regressions (with some bells and whistles) and the weights are automatically positive. That is, we do not need NMF at all! Once we have clusters and weights, we can study reconstruction accuracy and withincluster correlations between the underlying data and the fitted data that the cluster model produces.
We find that clustering works well for 10 out the 14 cancer types we study. The cancer types for which clustering does not appear to work all that well are Liver Cancer, Lung Cancer, and Renal Cell Carcinoma. Also, above 80% within-cluster correlations arise for 5 out of 7 clusters. Furthermore, remarkably, one cluster has high within-cluster correlations for 9 cancer types, and another cluster for 6 cancer types. These appear to be the leading clusters. Together they have high withincluster correlations in 11 out of 14 cancer types. So what does all this mean?
Additional insight is provided by looking at the within-cluster correlations between signatures Sig1 through Sig7 extracted in [8] and our clusters. High within-cluster correlations arise for Sig1, Sig2, Sig4 and Sig7, which are precisely the signatures with "peaks" (or "spikes" -"tall mountain landscapes"), whereas Sig3, Sig5 and Sig6 do not have such "peaks" ("flat" or "rolling hills landscapes"); see Figs. 14 through 20 of [8]. The latter 3 signatures simply do not have cluster-like structures. Looking at Fig. 21  In Section 2 we (i) discuss why applying clustering algorithms to extracting cancer signatures makes sense, (ii) argue that NMF, to a degree, is "clustering in disguise", and (iii) give the machinery for building cluster models via *K-means, including various details such as what to cluster, how to fix the number of clusters, etc. In Section 3 we discuss (i) cancer genome data we use, (ii) our application of *K-means to it, and (iii) the interpretation of our empirical results. Section 4 contains some concluding remarks, including a discussion of potential applications of *K-means in quantitative finance, where we outline some concrete problems where *K-means can be useful. Appendix A contains R source code for *K-means and cluster models.

Cluster models
The chief objective of this paper is to introduce a novel approach to identifying cancer signatures using clustering methods. In fact, as we discuss below in detail, our approach is more than just clustering. Indeed, it is evident from the get-go that blindly using nondeterministic clustering algorithms, 9 which typically produce (unmanageably) large numbers of local optima, would introduce great variability into the resultant cancer signatures. 10 On the other hand, deterministic algorithms such as agglomerative hierarchical clustering 11 typically are (substantially) slower and require essentially "guessing" the initial clustering, 12 which in practical applications 13 can often turn out to be suboptimal. So, both to motivate and explain our new approach employing clustering methods, we firstso to speak -"break down" the NMF approach and argue that it is in fact a clustering method in disguise!

"Breaking down" NMF
The current "lore"the commonly accepted method for extracting K cancer signatures from the occurrence counts matrix G is (see above) [5] is via nonnegative matrix factorization (NMF) [6,7]. Under NMF the matrix G is approximated via G ≈ W H, where W iA is an N × K matrix of weights, H As is a K × d matrix of exposures, and both W and H are nonnegative. However, not only is the number of signatures K not fixed via NMF (and must be either guessed or obtained via trial and error), NMF too is a nondeterministic algorithm and typically produces a large number of local optima. So, in practice one has no choice but to execute a large number N S of NMF runswhich we refer to as samplingsand then somehow extract cancer signatures from these samplings. Absent a guess for what K should be, one executes N S samplings for a range of values of K (say, K min ≤ K ≤ K max , where K min and K max are basically guessed based on some reasonable intuitive considerations), for each K extracts cancer signatures (see below), and then picks K and the corresponding signatures with the best overall fit into the underlying matrix G. For a given K, different samplings generally produce different weights matrices W. So, to extract a single matrix W for each value of K one averages over the samplings. However, before averaging, one must match the K cancer signatures across different samplingsindeed, in a given sampling X the columns in the matrix W iA are not necessarily aligned with the columns in the matrix W iA in a different sampling Y. To align the columns in the matrices W across the N S samplings, once often uses a clustering algorithm such as k-means. However, since k-means is nondeterministic, such alignment of the W columns is not guaranteed toand in fact does notproduce a unique answer. Here one can try to run multiple samplings of k-means for this alignment and aggregate them, albeit such aggregation itself would require another level of alignment (with its own nondeterministic clustering such as k-means). 14 And one can do this ad infinitum. In practice, one must break the chain at some level of alignment, either ad hoc (essentially by heuristically observing sufficient stability and "convergence") or via using a deterministic algorithm (see footnote 14 ). Either way, invariably all this introduces (overtly or covertly) systematic and statistical errors into the resultant cancer signatures and often it is unclear if they are meaningful without invoking some kind empirical biologic "experience" or "intuition" (often based on already well-known effects of, e.g., exposure to various well-understood carcinogens such as tobacco, ultraviolet radiation, aflatoxin, etc.). At the end of the day it all boils down to how usefulor predictivethe resultant method of extracting cancer signatures is, including signature stability. With NMF, the answer is not at all evident…

Clustering in disguise?
So, in practice, under the hood, NMF already uses clustering methods. However, it goes deeper than that. While NMF generically does not produce vanishing weights for a given signature, some weights are (much) smaller than others. E.g., often one has several "peaks" with high concentration of weights, with the rest of the mutation categories having relatively low weights. In fact, many weights can even be within the (statistical plus systematic) error bars. 15 Such weights can for all practical purposes be set to zero. In fact, we can take this further and ask whether proliferation of low weights adds any explanatory power. One way to address this is to run NMF with an additional constraint that the weights (obtained via averagingsee above) should be higher than either (i) some multiple of the corresponding error bars 16 or (ii) some preset fixed minimum weight. This certainly sounds reasonable, so why is this not done in practice? A prosaic answer appears to be that this would complicate the already nontrivial NMF algorithm even further, require additional coding and computation resources, etc. However, arguendo, let us assume that we require, say, that the weights be higher than a preset fixed minimum weight w min or else the weights are set to zero. As we increase w min , the so-modified NMF would produce more and more zeros. This does not mean that the resulting matrix W iA would have a binary cluster structure, i.e., that = W w δ i G i A iA ( ), , where δ AB is a Kronecker delta and G : {1, …, N} ↦ {1, …, K} is a map from N = 96 mutation categories to K clusters. Put another way, this does not mean that in the resulting matrix W iA for a given i (i.e., mutation category) we would have a nonzero element for one and only one value of A (i.e., signature). However, as we gradually increase w min , generally the matrix W iA is expected to look more and more like having a binary cluster structure, albeit with some "overlapping" signatures (i.e., such that in a given pair of signatures there are nonzero weights for one or more mutations). We can achieve a binary structure via a number of ways. Thus, a rudimentary algorithm would be to take the matrix W iA (equally successfully before or after achieving some zeros in it via nonzero w min ) and for a given value of i set all weights W iA to zero except in the signature A for which W iA = max(W iA |A = 1, …, K). Note that this might result in some empty signatures (clusters), i.e., signatures with W iA = 0 for all values of i. This can be dealt with by (i) ether simply dropping such signatures altogether and having fewer K′ < K signatures (binary clusters) at the end, or (ii) augmenting the algorithm to avoid empty clusters, which can be done in a number of ways we will not delve into here. The bottom line is that NMF essentially can be made into a clustering algorithm by reasonably modifying it, including via getting rid of ubiquitous and not-too-informative low weights. However, the downside would be an even more contrived algorithm, so this is not what we are suggesting here. Instead, we are observing that clustering is already intertwined in NMF and the question is whether we can simplify things by employing clustering methods directly.

Making clustering work
Happily, the answer is yes. Not only can we have much simpler and apparently more stable clustering algorithms, but they are also computationally much less costly than NMF. As mentioned above, the biggest issue with using popular nondeterministic clustering algorithms such as k-means 17 is that they produce a large number of local optima. For definiteness in the remainder of this paper we will focus on kmeans, albeit the methods described herein are general and can be applied to other such algorithms. Fortunately, this very issue has already been addressed in [11] in the context of constructing statistical industry classifications (i.e., clustering models for stocks) for quantitative trading, so here we simply borrow therefrom and further expand and adapt that approach to cancer signatures.

K-means
A popular clustering algorithm is k-means [12][13][14][15][16][17][18]. The basic idea behind k-means is to partition N observations into K clusters such that each observation belongs to the cluster with the nearest mean. Each of the N observations is actually a d-vector, so we have an N × d matrix X is , i = 1, …, N, s = 1, …, d. Let C a be the K clusters, C a = {i|i ∈ C a }, a = 1, …, K. Then k-means attempts to minimize 18 are the cluster centers (i.e., cross-sectional means), 19 and n a = |C a | is the number of elements in the cluster C a . In (1) the measure of "closeness" is chosen to be the Euclidean distance between points in R d , albeit other measures are possible.
One "drawback" of k-means is that it is not a deterministic algorithm. Generically, there are copious local minima of g in (1) and the algorithm only guarantees that it will converge to a local minimum, not the global one. Being an iterative algorithm, unless the initial centers are preset, k-means starts with a random set of the centers Y as at the initial iteration and converges to a different local minimum in each run. There is no magic bullet here: in practical applications, typically, trying to "guess" the initial centers is not any easier than "guessing" where, e.g., the global minimum is. So, what is one to do? One possibility is to simply live with the fact that every run produces a different answer. In fact, this is acceptable in many applications. However, in the context of extracting cancer signatures this would result in an exercise in futility. We need a way to eliminate or greatly reduce indeterminism.

Aggregating clusterings
The idea is simple. What if we aggregate different clusterings from multiple runswhich we refer to as samplingsinto one? The question is how. Suppose we have M runs (M ≫ 1). Each run produces a clustering with K clusters. Let ). Here we are assuming that somehow we know how to properly order (i.e., align) the K clusters from each run. This is a nontrivial assumption, which we will come back to momentarily. However, assuming, for a second, that we know how to do this, we can aggregate the binary matrices Ω r ia into a single matrix ia . Now, this matrix does not look like a binary clustering matrix. Instead, it is a matrix of occurrence counts, i.e., it counts how many times a given mutation was assigned to a given cluster in the process of M samplings. What we need to construct is a map G such that one and only one mutation belongs to each of the K clusters. The simplest criterion is to map a given mutation to the cluster in which ∼ Ω ia is maximal, i.e., where said mutation occurs most frequently. A caveat is that there may be more than one such clusters. A simple criterion to resolve such an ambiguity is to assign said mutation to the cluster with most cumulative occurrences (i.e., we assign said mutation to the cluster with the largest ia ). Further, in the unlikely event that there is still an ambiguity, we can try to do more complicated things, or we can simply assign such a mutation to the cluster with the lowest value of the index atypically, there is so much noise in the system that dwelling on such minutiae simply does not pay off.
However, we still need to tie up a loose end, to wit, our assumption that the clusters from different runs were somehow all aligned. In practice each run produces K clusters, but (i) they are not the same clusters and there is no foolproof way of mapping them, especially when we have a large number of runs; and (ii) even if the clusters were the same or similar, they would not be ordered, i.e., the clusters from one run generally would be in a different order than the clusters from another run.
So, we need a way to "match" clusters from different samplings. Again, there is no magic bullet here either. We can do a lot of complicated and contrived things with not much to show for it at the end. A simple pragmatic solution is to use k-means to align the clusters from different runs. Each run labeled by r = 1, …, M, among other things, produces a set of cluster centers Y r as . We can "bootstrap" them by row into a (KM) takes values 17 Which are preferred over deterministic ones for the reasons discussed above. 18 Below we will discuss what X is should be for cancer signatures. 19 Throughout this paper "cross-sectional" refers to "over the index i". 20 Note that here the superscript r in Ω r ia , G r (i) and N a r (see below) is an index, not a power.
Z. Kakushadze, W. Yu Biomolecular Detection and Quantification 13 (2017) 7-31 = … a 1, , (KM) ͠ . We can now cluster ∼ Y as ͠ into K clusters via k-means. This will map each value of a͠ to {1, …, K} thereby mapping the K clusters from each of the M runs to {1, …, K}. So, this way we can align all clusters. The "catch" is that there is no guarantee that each of the K clusters from each of the M runs will be uniquely mapped to one value in {1, …, K}, i.e., we may have some empty clusters at the end of the day. However, this is fine, we can simply drop such empty clusters and aggregate (via the above procedure) the smaller number of K′ < K clusters. I.e., at the end we will end up with a clustering with K′ clusters, which might be fewer than the target number of clusters K. This is not necessarily a bad thing. The dropped clusters might have been redundant in the first place. Another evident "catch" is that even the number of resulting clusters K′ is not deterministic. If we run this algorithm multiple times, we will get varying values of K′. Malicious circle?

Fixing the "ultimate" clustering
Not really! There is one other trick up our sleeves we can use to fix the "ultimate" clustering thereby rendering our approach essentially deterministic. The idea above is to aggregate a large enough number M of samplings. Each aggregation produces a clustering with some K′ ≤ K clusters, and this K′ varies from aggregation to aggregation. However, what if we take a large number P of aggregations (each based on M samplings)? Typically there will be a relatively large number of different clusterings we get this way. However, assuming some degree of stability in the data, this number is much smaller than the number of a priori different local minima we would obtain by running the vanilla kmeans algorithm. What is even better, the occurrence counts of aggregated clusterings are not uniform but typically have a (sharply) peaked distribution around a few (or manageable) number of aggregated clusterings. In fact, as we will see below, in our empirical genome data we are able to pinpoint the "ultimate" clustering! So, to recap, what we have done here is this. There are myriad clusterings we can get via vanilla k-means with little to no guidance as to which one to pick. 21 We have reduced this proliferation by aggregating a large number of such clusterings into our aggregated clusterings. We then further zoom onto a few or even a unique clustering we consider to be the likely "ultimate" clustering by examining the occurrence counts of such aggregated clusterings, which turns out to have a (sharply) peaked distribution. Since vanilla k-means is a relatively fast-converging algorithm, each aggregation is not computationally taxing and running a large number of aggregations is nowhere as time consuming as running a similar number (or even a fraction thereof) of NMF computations (see below).

What to cluster?
So, now that we know how to make clustering work, we need to decide what to cluster, i.e., what to take as our matrix X is in (1). The naïve choice X is = G is is suboptimal for multiple reasons (as discussed in [8]).
First, the elements of the matrix G is are populated by nonnegative occurrence counts. Nonnegative quantities with large numbers of samples tend to have skewed distributions with long tails at higher values. I.e., such distributions are not normal but (in many cases) roughly log-normal. One simple way to deal with this is to identify X is with a (natural) logarithm of G is (instead of G is itself). A minor hiccup here is that some elements of G is can be 0. We can do a lot of complicated and even convoluted things to deal with this issue. Here, as in [8], we will follow a pragmatic approach and do something simple insteadthere is so much noise in the data that doing convoluted things simply does not pay off. So, as the first cut, we can take This takes care of the G is = 0 cases; for G is ≫ 1 we have R is ≈ ln(G is ), as desired. Second, the detailed empirical analysis of [8] uncovered what is termed therein the "overall" mode 22 unequivocally present in the occurrence count data. This "overall" mode is interpreted as somatic mutational noise unrelated to (and in fact obscuring) the true underlying cancer signatures and must therefore be factored out somehow.
Here is a simple way to understand the "overall" mode. Let the correlation matrix Ψ ij = Cor(X is , X js ), where Cor(·, ·) is serial correlation. 23 I.e., ii are variances, and the serial covariance matrix 24 ij between different mutation categories is nonzero and is in fact high for most cancer types we study. This is the aforementioned somatic mutational noise that must be factored out. If we aggregate samples by cancer types (see below) and compute the correlation matrix Ψ ij for the so-aggregated data (across the n = 14 cancer types we studysee below), 25 the average correlation ρ is over whopping 96%. Another way of thinking about this is that the occurrence counts in different samples (or cancer types, if we aggregate samples by cancer types) are not normalized uniformly across all samples (cancer types). Therefore, running NMF, a clustering or any other signature-extraction algorithm on the vanilla matrix G is (or its "log" X is defined in (3)) would amount to mixing apples with oranges thereby obscuring the true underlying cancer signatures.
Following [8], factoring out the "overall" mode (or "de-noising" the matrix G is ) therefore most simply amount to cross-sectional (i.e., across the 96 mutation categories) demeaning of the matrix X is . I.e., instead of X is we use ′ X is , which is obtained from X is by demeaning its columns: 26 We should note that using ′ X is instead of X is in (1) does not affect clustering. Indeed, g in (1) is invariant under the transformations of the form X is → X is + Δ s , where Δ s is an arbitrary d-vector, as thereunder we also have Y as → Y as + Δ s , so X is − Y as is unchanged. In fact, this is good: this means that de-noising does not introduce any additional errors into clustering itself. However, the actual weights in the matrix W iA are affected by de-noising. We discuss the algorithm for fixing W iA below. However, we need one more ingredient before we get to determining the weights, and with this additional ingredient de-noising does affect clustering.

Normalizing log-counts
As was discussed in [11], clustering X is (or equivalently ′ X is ) would be suboptimal. 27 The issue is this. Let ′ σ i be serial standard deviations, references therein) corresponding to the overall movement of the broad market, which affects all stocks (to varying degrees)cash inflow (outflow) into (from) the market tends to push stock prices higher (lower). This is the market risk factor, and to mitigate it one can, e.g., hold a dollar-neutral portfolio of stocks (i.e., the same dollar holdings for long and short positions). 23 Throughout this paper "serial" refers to "over the index s". 24 The overall normalization of C ij , i.e., d − 1 (unbiased estimate) vs. d (maximum likelihood estimate) in the denominator in the definition of C ij in (4), is immaterial for our purposes here. 25 So, in this case d = n = 14 in (4). 26 For the reasons discussed above, we should demean X is , not G is . 27 More precisely, the discussion of [11] is in the financial context, to wit, quantitative trading, which has its own nuances (see below). However, some of that discussion is quite is , where, as above, Cov(·, ·) is serial covariance. Here we assume that samples are aggregated by cancer types, so s = 1, …, d with d = n = 14. Now, ′ σ i are not cross-sectionally uniform and vary substantially across mutation categories. The density of ′ σ i is depicted in Fig. 1  X is , this variability in ′ σ i will not be accounted for.
A simple solution is to cluster normalized demeaned log-counts This way we factor out the nonuniform (and skewed) standard deviation out of the log-counts. Note that now denoising does make a difference in clustering. Indeed, if we use (1) and (2), the quantity g (and also clusterings) will be different.

Fixing cluster number
Now that we know what to cluster (to wit, ′ ∼ X is ) and how to get to the "unique" clustering, we need to figure out how to fix the (target) number of clusters K, which is one of the inputs in our algorithm above. 29 In [8] it was argued that in the context of cancer signatures their number can be fixed by building a statistical factor model [9], i.e., the number of signatures is simply the number of statistical factors. 30 So, by the same token, here we identify the (target) number of clusters in our clustering algorithm with the number of statistical factors fixed via the method of [9].

Effective rank
So, following [9,8], we set 31 Here eRank(Z) is the effective rank [10] of a symmetric semi-positivedefinite (which suffices for our purposes here) matrix Z. It is defined as (8)(9) where λ (a) are the L positive eigenvalues of Z, and H has the meaning of the (Shannon a.k.a. spectral) entropy [34,35]. Let us emphasize that in (6) the matrix Ψ ij is computed based on the demeaned log-counts 32 ′ X is . The meaning of eRank(Ψ ij ) is that it is a measure of the effective dimensionality of the matrix Ψ ij , which is not necessarily the same as the number L of its positive eigenvalues, but often is lower. This is due to the fact that many d-vectors ′ X is can be serially highly correlated (which manifests itself by a large gap in the eigenvalues) thereby further reducing the effective dimensionality of the correlation matrix.

How to compute weights?
The one remaining thing to accomplish is to figure out how to compute the weights W iA . Happily, in the context of clustering we have significant simplifications compared with NMF and computing the weights becomes remarkably simple once we fix the clustering, i.e., the matrix Ω iA = δ G(i),A (or, equivalently, the map G : {i} ↦ {A}, i = 1, …, N, A = 1, …, K, where for the notational convenience we use K to denote the number of clusters in the "ultimate" clusteringsee above). Just as in NMF, we wish to approximate the matrix G is via a product of the weights matrix W iA and the exposure matrix H As , both of which must be nonnegative. More precisely, since we must remove the "overall" mode, i.e., de-noise the matrix G is , following [8], instead of G is we will approximate the re-exponentiated demeaned log-counts matrix ′ X is : We can include an overall normalization by taking s is is (recall that X s is the vector of column means of X issee Eq. (5)), etc., to make it look more like the original matrix G is ; however, this does not affect the extracted signatures. 33 Also, technically speaking, after re-exponentiating we should "subtract" the extra 1 we added in the definition (3) (assuming we include one of the aforesaid overall normalizations). However, the inherent noise in the data makes this a moot point.
So, we wish to approximate ′ G is via a product W H. However, with clustering we have = W w δ i G i A iA ( ), , i.e., we have a block (cluster) structure where for a given value of A all W iA are zero except for i ∈ J Fig. 1. Horizontal axis: serial standard deviation ′ σ i for N = 96 mutation categories (i = 1, …, N) of cross-sectionally demeaned log-counts ′ X is across n = 14 cancer types (for samples aggregated by cancer types, so s = 1, …, d, d = n). Vertical axis: density using R function density(). See Section 2.4.1 for details.
(footnote continued) general and can be adapted to a wide variety of applications. 28 Qu. = Quartile, SD = Standard Deviation, MAD = Mean Absolute Deviation. 29 A variety of methods for fixing the number of clusters have been discussed in other contexts, e.g., [22][23][24][25][26][27][28][29]. 30 In the financial context, these are known as statistical risk models [9]. For a discussion and literature on multifactor risk models, see, e.g., [30,31] and references therein. For prior works on fixing the number of statistical risk factors, see, e.g., [32,33]. 31 Here Round(·) can be replaced by floor(·) = ⌊·⌋. 32 Note that using normalized demeaned log-counts ′ ∼ X is gives the same Ψ ij . 33 This is because each column of W, being weights, is normalized to add up to 1.

Z. Kakushadze, W. Yu
Biomolecular Detection and Quantification 13 (2017)  (A) = {j|G(j) = A}, i.e., for the mutation categories labeled by i that belong to the cluster labeled by A. Therefore, our matrix factorization of G is into a product W H now simplifies into a set of K independent factorizations as follows: So, there is no need to run NMF anymore! Indeed, if we can somehow fix H As for a given cluster, then within this cluster we can determine the corresponding weights w i (i ∈ J(A)) via a serial linear regression: where ε is are the regression residuals. I.e., for each A ∈ {1, …, K}, we 2 . Furthermore, if H As are nonnegative, then the weights w i are automatically nonnegative as they are given by: Now, we wish these weights to be normalized: This can always be achieved by rescaling H As . Alternatively, we can pick H As without worrying about the normalization, compute w i via (13), rescale them so that they satisfy (14), and simultaneously accordingly rescale H As . Mission accomplished!

Fixing exposures
Well, almost… We still need to figure out how to fix the exposures H As . The simplest way to do this is to note that we can use the matrix Ω iA = δ G(i),A to swap the index i in ′ G is by the index A, i.e., we can take That is, up to the normalization constants η͠ A (which are fixed via (14)) we simply take cross-sectional means of ′ G is in each cluster. (Recall that n A = J(A).) The so-defined H As are automatically positive as all ′ G is are positive. Therefore, w i defined via (13) are also all positive. This is a good newsvanishing w i would amount to an incomplete weights matrix W iA (i.e., some mutations would belong to no cluster).
So, why does (15) make sense? Looking at (12), we can observe that, if the residuals ε is cross-sectionally, within each cluster labeled by A, are random, then we expect that ∑ i∈J(A) ε is ≈ 0. If we had an exact equality here, then we would have (15) with η A = 1 (i.e., = η n ͠ A A ) assuming the normalization (14). In practice, the residuals ε is are not exactly "random". First, the number n A of mutation categories in each cluster is not large. Second, as mentioned above, there is variability in serial standard deviations across mutation types. This leads us to consider variations.

A variation
Above we argued that it makes sense to cluster normalized demeaned log-counts ′ = ′ ′ ∼ X X σ / i is is due to the cross-sectional variability (and skewness) in the serial standard deviations ′ σ i . We may worry about similar effects in ′ G is when computing H As and w i as we did above. This can be mitigated by using normalized quantities is are serial variances. That is, we can define 35 where ν A = ∑ i∈J(A) 1/ω i . So, 1/ω i are the weights in the averages over the clusters.

Another variation
Here one may wonder, considering the skewed roughly log-normal distribution of G is and henceforth ′ G is , would it make sense to relate the exposures to within-cluster cross-sectional averages of demeaned logcounts ′ X is as opposed to those of ′ G is ? This is easily achieved. Thus, we can define (this ensures positivity of H As ): Exponentiating we get I.e., instead of an arithmetic average as in (15), here we have a geometric average. As above, here too we can introduce nontrivial weights. Note that the form of (17) is the same as (13), it is only H As that is affected by the weights. So, we can introduce the weights in the geometric means as follows: is . Thus, we have: So, the weights are the exponents ′ μ σ 1/ A i . Other variations are also possible.

Implementation
We are now ready to discuss an actual implementation of the above algorithm, much of the R code for which is already provided in [8,11]. The R source code is given in Appendix A hereof.

Data summary
In our empirical analysis below we use the same genome data (from published samples only) as in [8]. This data is summarized in Table S1 (borrowed from [8]), which gives total counts, number of samples and the data sources, which are as follows: A1 = [36] [52]. Sample IDs with the corresponding publication sources are given in Appendix A of [8]. In our analysis below we aggregate samples by the 14 cancer types. The resulting data is in Tables S2 and S3. For tables and figures labeled S★ see Supplementary Materials (see Appendix C for a web link).

Structure of data
The underlying data consists of a matrixcall it G iswhose elements are occurrence counts of mutation types labeled by i = 1, …, N = 96 in samples labeled by s = 1, …, d. More precisely, we can work 34 The superscript T denotes matrix transposition. 35 I.e., here we assume that ε is /ω i are approximately random in (12). with one matrix G is which combines data from different cancer types; or, alternatively, we may choose to work with individual matrices [G(α)] is , where: α = 1, …, n labels n different cancer types; as before, i = 1, …, N = 96; and s = 1, …, d(α). Here d(α) is the number of samples for the cancer type labeled by α. The combined matrix G is is obtained simply by appending (i.e., bootstrapping) the matrices [G(α)] is together column-wise. In the case of the data we use here (see above), this "big matrix" turns out to have 1389 columns. Generally, individual matrices [G(α)] is and, thereby, the "big matrix", contain a lot of noise. For some cancer types we can have a relatively small number of samples. We can also have "sparsely populated" data, i.e., with many zeros for some mutation categories. As mentioned above, different samples are not necessarily uniformly normalized. Etc. The bottom line is that the data is noisy. Furthermore, intuitively it is clear that the larger the matrix we work with, statistically the more "signatures" (or clusters) we should expect to get with any reasonable algorithm. However, as mentioned above, a large number of signatures would be essentially useless and defy the whole purpose of extracting them in the first placewe have 96 mutation categories, so it is clear that the number of signatures cannot be more than 96! If we end up with, say, 50+ signatures, what new or useful does this tell us about the underlying cancers? The answer is likely nothing other than that most cancers have not much in common with each other, which would be a disappointing result from the perspective of therapeutic applications. To mitigate the aforementioned issues, at least to a certain extent, following [8], we can aggregate samples by cancer types. This way we get an N × n matrix, which we also refer to as G is , where the index s = 1, …, d now takes d = n values corresponding to the cancer types. In the data we use n = 14, the aggregated matrix G is is much less noisy than the "big matrix", and we are ready to apply the above machinery to it.

Genome data results
The 96 × 14 matrix G is given in Tables S2 and S3 is what we pass into the function bio.cl.sigs() in Appendix A as the input matrix x. We use: iter.max = 100 (this is the maximum number of iterations used in the built-in R function kmeans()we note that there was not a single instance in our 150 million runs of kmeans() where more iterations were required); 36 num.try = 1000 (this is the number of individual k-means samplings we aggregate every time); and num.runs = 150000 (which is the number of aggregated clusterings we use to determine the "ultimate"that is, the most frequently occurringclustering). So, we ran k-means 150 million times. More precisely, we ran 15 batches with num.runs = 10000 as a sanity check, to make sure that the final result based on 150,000 aggregated clusterings was consistent with the results based on smaller batches, i.e., that it was insample stable. 37 Based on Table S4, we identify Clustering-A as the "ultimate" clustering (cf. Clustering-B/C/D).
We give the weights for Clustering-A, Clustering-B, Clustering-C and Clustering-D using unnormalized and normalized regressions with exposures computed based on arithmetic averages (see Section 2.6) in Tables 1, 2, S5-S10, and Figs. 2 through Fig. 15 and S1 through S40. We give the weights for Clustering-A using unnormalized and normalized regressions with exposures computed based on geometric averages (see Section 2.6) in Tables 3, 4, and Figs. S41 through S54. The actual mutation categories in each cluster for a given clustering can be read off the aforesaid tables with the weights (the mutation categories with nonzero weights belong to a given cluster), or from the horizontal axis labels in the aforesaid figures. It is evident that Clustering-A, Clustering-B, Clustering-C and Clustering-D are essentially variations of each other (Clustering-D has only 6 clusters, while the other 3 have 7 clusters).

Reconstruction and correlations
So, based on genome data, we have constructed clusterings and weights. Do they work? I.e., do they reconstruct the input data well? It is evident from the get-go that the answer to this question may not be binary in the sense that for some cancer types we might have a nice clustering structure, while for others we may not. The aim of the following exercise is to sort this all out. Here come the correlations…

Within-cluster correlations
We have our de-noised 38 matrix ′ G is . We are approximating this matrix via the following factorized matrix: We can now compute an n × K matrix Θ sA of within-cluster cross-sectional correlations between ′ G is and G* is defined via (xCor(·, ·) stands for "cross-sectional correlation" to distinguish it from "serial correlation" Cor(·, ·) we use above) 39 We give this matrix for Clustering-A with weights using normalized regressions with exposures computed based on arithmetic means (see Section 2.6) in Table 5. Let us mention that, with exposures based on arithmetic means, weights using normalized regressions work a bit better than using unnormalized regressions. Using exposures based on geometric means changes the weights a bit, which in turn slightly affects the within-cluster correlations, but does not alter the qualitative picture.

Overall correlations
Another useful metric, which we use as a sanity check, is this. For each value of s (i.e., for each cancer type), we can run a linear crosssectional regression (without the intercept) of ′ G is over the matrix W iA . So, we have n = 14 of these regressions. Each regression produces multiple R 2 and adjusted R 2 , which we give in Table 5. Furthermore, we can compute the fitted valuesĜ * is based on these regressions, which are given by where (for each value of s) F As are the regression coefficients. We can now compute the overall cross-sectional correlations (i.e., the index i runs over all N = 96 mutation categories) These correlations are also given in Table 5 and measure the overall fit quality.

Interpretation
Looking at Table 5 a few things become immediately evident. Clustering works well for 10 out the 14 cancer types we study here. The cancer types for which clustering does not appear to work all that well are Breast Cancer (labeled by X4 in Table 5), Liver Cancer (X8), Lung 36 The R function kmeans() produces a warning if it does not converge within iter.max. 37 We ran these 15 batches consecutively, and each batch produced the same top-10 (by occurrence counts) clusterings as in Table S4; however, the actual occurrence counts are different across the batches with slight variability in the corresponding rankings. The results are pleasantly stable. 38 De-noising per se does not affect cross-sectional correlations. Adding extra 1 in (3) (recall that we obtain ′ G is by cross-sectionally demeaning X is and then re-exponentiating) has a negligible effect. So, in the correlations below we can use the original data matrix G is instead of ′ G is . 39 Due to the factorized structure (22), these correlations do not directly depend on H As .

Z. Kakushadze, W. Yu
Biomolecular Detection and Quantification 13 (2017)  Cancer (X9), and Renal Cell Carcinoma (X14). More precisely, for Breast Cancer we do have a high within-cluster correlation for Cl-5 (and also Cl-4), but the overall fit is not spectacular due to low within-cluster correlations in other clusters. Also, above 80% within-cluster correlations 40 arise for 5 clusters, to wit, Cl-1, Cl-3, Cl-4, Cl-5 and Cl-6, but not for Cl-2 or Cl-7. Furthermore, remarkably, Cl-1 has high within-cluster correlations for 9 cancer types, and Cl-5 for 6 cancer types. These appear to be the leading clusters. Together they have high within-cluster correlations in 11 cancer types. So what does all this mean? Additional insight is provided by looking at the within-cluster correlations between the 7 cancer signature extracted in [8] and the clusters we find here. Let iα W be the weights for the 7 cancer signatures from Tables 13 and 14 of [8]. We can compute the following withincluster correlations (α = 1, …, 7 labels the cancer signatures of [8], which we refer to as Sig1 through Sig7): These correlations are given in Table 6. High within-cluster correlations arise for Cl-1 (with Sig1 and Sig7), Cl-5 (with Sig2) and Cl-6 (with Sig4). And this makes perfect sense. Indeed, looking at Figs. 14 through 20 of [8], Sig1, Sig2, Sig4 and Sig7 are precisely the cancer signatures that have "peaks" (or "spikes" -"tall mountain landscapes"), whereas Sig3, Sig5 and Sig6 do not have such "peaks" ("flat" or "rolling hills landscapes"). No wonder such signatures do not have high withincluster correlationsthey simply do not have cluster-like structures. Looking at Fig. 21 in [8], it becomes evident why clustering does not work well for Liver Cancer (X8)it has a whopping 96% contribution from Sig5! Similarly, Renal Cell Carcinoma (X14) has a 70% contribution from Sig6. Lung Cancer (X9) is dominated by Sig3, hence no cluster-like structure. Finally, Breast Cancer (X4) is dominated by Sig2,  Table 2 with the next 48 mutation categories) for the 7 clusters in Clustering-A (see Table S4) based on unnormalized (columns 2-8) and normalized (columns 9-15) regressions (see Section 2.6 for details). Each cluster is defined as containing the mutations with nonzero weights. (The mutations are encoded as follows: XYZW = Y>W: XYZ. Thus, GCGA = C>A: GCG.) For instance, cluster Cl-2 contains 8 mutations GCGA, TCGA, ACGG, GCCG, GCGG, TCGG, GTCA, GTCG. In each cluster the weights are normalized to add up to 100% (up to 2 digits due to the aforesaid rounding). In Tables 1 through S10 "weights based on unnormalized regressions" are given by (13), (14) and (15), while "weights based on normalized regressions" are given by (17), (14) and (16), i.e., the exposures are calculated based on arithmetic averages (see Section 2.6 for details). which has a high within-cluster correlation with Cl-5, which is why Breast Cancer has a high within-cluster correlation with Cl-5 (but poor overall correlation in Table 5 Cancer (X13) apparently all have the Cl-5 cluster structure embedded in them substantially. Furthermore, note the overlap between these two lists, to wit, Esophageal Cancer(X6), Ovarian Cancer (X11), Pancreatic Cancer (X12) and Prostate Cancer (X13). We obtained this result purely statistically, with no biologic input, using our clustering algorithm and other statistical methods such as linear regression to obtain the actual weights. It is too early to know whether this insight will aid any therapeutic applications, but that is the hopesimilarities in the underlying genomic structures of different cancer types raise hope that therapeutics for one cancer type could perhaps be applicable to other cancer types. On the other hand, our findings above relating to Liver Cancer, Lung Cancer and Renal Cell Carcinoma (and possibly also Breast Cancer, albeit the latter does appear to have a not-so-insignificant overlap with Cl-5, which differentiates it from the aforesaid 3 cancer types) suggest that these cancer types apparently stand out.

Concluding remarks
Clustering ideas and techniques have been applied in cancer research in various incarnations and contexts aplentyfor a partial list of works at least to some extent related to our discussion here, see, e.g., [52,53,54,55,40,56,5,36, and references therein. As mentioned above, even in NMF clustering is used at some (perhaps not-so-evident) layer. What is new in our approachand hence new resultsis that: (i) following [8], we apply clustering to aggregated by cancer types and   de-noised data; ii) we use a tried-and-tested in quantitative finance bag of tricks from [11], which improves clustering; and (iii) last but not least, we apply our *K-means algorithm to cancer genome data. As mentioned above, *K-means, unlike vanilla k-means or its other commonly used variations, is essentially deterministic, and it achieves determinism statistically, not by "guessing" initial centers or as in  agglomerative hierarchical clustering, which basically "guesses" the initial (e.g., 2-cluster) clustering. Instead, via aggregating a large number of k-means clusterings and statistical examination of the occurrence counts of such aggregations, *K-means takes a mess of myriad vanilla k-means clusterings and systematically reduces randomness and indeterminism without ad hoc initial "guesswork".      As mentioned above, consistently with the results of [8] obtained via improved NMF techniques, Liver Cancer, Lung Cancer and Renal Cell Carcinoma do not appear to have clustering (sub)structures. This could be both good and bad news. It is a good news because we learned something interesting about these cancer typesand in two complementary ways. However, it could also be a bad news from the therapeutic standpoint. Since these cancer types appear to have little in common with others, it is likely that they would require specialized therapeutics. On the flipside, we should note that it would make sense to exclude these 3 cancer types when running clustering analysis. However, it would also make sense to include other cancer types by utilizing the International Cancer Genome Consortium data, which we leave for future studies. (For comparative reasons, here we used the same data as in [8], which was limited to data samples published as of the date thereof.) This paper is not intended to be an exhaustive empirical study but a proof of concept and an opening of a new avenue for extracting and studying cancer signatures beyond the tools that NMF provides.
And we do find that 11 out of the 14 cancer types we study here have clustering structures substantially embedded in them and clustering overall works well for at least 10 out of these 11 cancer types. 41 Now, looking at Fig. 14 of [8], we see that its "peaks" are located at ACGT, CCGT, GCGT and TCGT. The same "peaks" are present in our cluster Cl-1 (see Figs. 2 and 3). Hence the high within-cluster correlation between Cl-1 and Sig1. On the other hand, Sig1 of [8] is essentially the same as the mutational signature 1 of [40,36], which is due to spontaneous cytosine deamination. So, this is what our cluster Cl-1 describes. Next, looking at Fig. 15 of [8], we see that its "peaks" are located at TCAG, TCTG, TCAT and TCTT. The first two of these "peaks" TCAG and TCTG are present in our Cl-5 (see Figs. 10 and 11), the third "peak" TCAT is present in our Cl-1 (see Figs. 2 and 3), while the fourth "peak" TCTT is present in our Cl-4 (see Figs. 8 and 9), which is consistent with the high within-cluster correlations between Sig2 and Cl-4 Table 3 Weights (in the units of 1%, rounded to 2 digits) for the first 48 mutation categories for the 7 clusters in Clustering-A (see Table S4) based on unnormalized (columns 2-8) and normalized (columns 9-15) regressions with the exposures computed via geometric means (see Section 2.6 for details). Here "weights based on unnormalized regressions" are given by (13), (14) and (19), while "weights based on normalized regressions" are given by (17), (14) and (21). Other conventions are the same as in Table 1.   Table 5 The within-cluster cross-sectional correlations Θ sA (columns 2-8), the overall correlations Ξ s (column 11) based on the overall cross-sectional regressions, and multiple R 2 and adjusted R 2 of these regressions (columns 9 and 10). See Section 3.3 for details. Cancer types are labeled by X1 through X14 as in and Cl-5, albeit its within-cluster correlation with Cl-1 is poor. Note that Sig2 of [8] is essentially the same as the mutational signatures 2 + 13 of [40,36], which are due to APOBEC mediated cytosine deamination. In fact, it was reported as a single signature in [36], however, subsequently, it was split into 2 distinct signatures, which usually appear in the same samples. 42 Our clustering results indicate that grouping TCAG and TCTG into one signature makes sense as they belong to the same cluster Cl-5. However, grouping TCAT and TCTT together does not appear to make much sense. Looking at the figures for Clustering-A, Clustering-B, Clustering-C and Clustering-D, we see that the TCAT "peak" invariably appears together with the ACGT, CCGT, GCGT and TCGT "peaks" as in Cl-1 in Clustering-A, Cl-2 in Clustering-B, Cl-1 in Clustering-C, and Cl-1 in Clustering-D, but never with TCTT. So, our clustering approach tells us something new beyond the NMF "intuition". This may have an important implication for Breast Cancer, which, as mentioned above, is dominated by Sig2. Thus, based on our results in Table 5, we see that Breast Cancer has high within-cluster correlations with Cl-4 and Cl-5, but not with Cl-1. This may imply that clustering simply does not work well for Breast Cancer, which would appear to put it in the same "stand-alone" league as Liver Cancer, Lung Cancer and Renal Cell Carcinoma. In any event, clustering invariably suggests that the TCAT "peak" belongs in Cl-1 with the 4 "peaks" ACGT, CCGT, GCGT and TCGT related to spontaneous cytosine deamination, rather than those related to APOBEC mediated cytosine deamination. Now, let us check the remaining two signatures of [8] with "tall mountain landscapes" (see above), to wit, Sig4 and Sig7. Looking at Fig. 17 of [8], we see that its "peaks" are at CTTC, TTTC, CTTG and TTTG. The same peaks appear in our Cl-6 (see Figs. 12 and 13). Hence the high within-cluster correlation between Cl-6 and Sig4. Note that Sig4 is essentially the same as the mutational signature 17 of [40,36], and its underlying mutational process is unknown. Next, looking at Fig. 20 of [8], we see that its "peaks" for the C > G mutations are essentially the same as in Cl-1. Hence the high within-cluster correlation between Cl-7 and Sig1. So, there are no surprises with Sig1, Sig4 and Sig7. However, based on our clustering results, as we discuss above, with Sig2 we do findwhat we feel is a pleasantsurprise, that splitting it into two signatures (see above) might be inadequate and the TCAT "peak" might really belong with the Sig1 "peaks" (spontaneous v. APOBEC mediated cytosine deamination). This is exciting as it might be an indication of the limitations of NMF (or clustering…). 43 In Introduction we promised that we would discuss some potential applications of *K-means in quantitative finance, and so here it is. Let us mention that *K-means is universal, oblivious to the input data and applicable in a variety of fields. In quantitative finance *K-means a priori can be applied everywhere clustering methods are used with the added bonus of (statistical) determinism. 44 One evident example is statistical industry classifications discussed in [11], where one uses clustering methods to classify stocks. In fact, *K-means is an extension of the methods discussed in [11]. One thing to keep in mind is that in *K-means one sifts through a large number P of aggregations, which can get computationally costly when clustering 2000+ stocks into 100+ clusters. 45 Another potential application is in the context of combining alphas (trading signals)see, e.g., [79]. Yet another application is when we have a term structure, such as a portfolio of bonds (e.g., U.S. Treasuries or some other bonds) with varying maturities, or futures (e.g., Eurodollar futures) with varying deliveries. These cases resemble the genome data more in the sense that the number N of instruments is relatively small (typically even fewer than the number of mutation categories). Another example with a relatively small number of instruments would be a portfolio of various futures for different FX (foreign exchange) pairs (even with the uniform delivery), e.g., USD/ EUR, USD/HKD, EUR/AUD, etc., i.e., FX statistical arbitrage. One approach to optimizing risk in such portfolios is by employing clustering methods and a stable, essentially deterministic algorithm such as *Kmeans can be useful. Hopefully *K-means will prove a valuable tool in cancer research, quantitative finance as well as various other fields (e.g., image recognition).

Conflict of interest
Authors declare no conflict of interest.

Appendix A. R source code
In this appendix we give the R (R Package for Statistical Computing, http://www.r-project.org) source code for computing the clusterings and weights using the algorithms of Section 2. The code is straightforward and self-explanatory. 46 The main function is bio.cl.sigs(x, iter.max = 100, num.try = 1000, num.runs = 10000). Here: x is the N × d occurrence counts matrix G is (where N = 96 is the number of mutation categories, and d is the number of samples; or d = n, where n is the number of cancer types, when the samples are aggregated by cancer types); iter.max is the maximum number of iterations that are passed into the R built-in function kmeans(); num.try is the number M of aggregated clusterings (see Section 2.3.2); num.runs is the number of runs P used to determine the most frequently occurring clustering (the "ultimate" clustering) obtained via aggregation (see Section 2.3.3). The function bio.erank.pc() is defined in Appendix B of [8]. The function qrm.stat.ind.class() is defined in Appendix A of [11]. This function internally calls another function qrm.calc.norm.ret(), which we redefine here via the function bio.calc.norm.ret(). 47 The output is a list, whose elements are as follows: res$ind is an N × K binary matrix Ω iA = δ G(i),A (i = 1, …, N, A = 1, …, K, the map G : {1, …, N} ↦ {1, …, K}see Section 2), which defines the K clusters in the "ultimate" clustering; 48 res$w is an N-vector of weights obtained via unnormalized regressions using arithmetic means for computing Table 6 The within-cluster cross-sectional correlations Δ αA between the weights for 7 cancer signatures Sig1 through Sig7 of [8] and the weights (using normalized regressions with exposures based on arithmetic averages) for 7 clusters in Clustering A (see Section 3.3 for details). All quantities are in the units of 1% rounded to 2 digits. The values above 80% are given in bold font.  43 Or both… Alternativelyand that would be truly excitingperhaps there is a biologic explanation. In any event, it is too early to tellyet another possibility is that this is merely an