Copyright for this paper by its authors. Use permitted under Creative Commons
License Attribution 4.0 International (CC BY 4.0).

Computational Humanities Research Conference 2021

[orcid=0000-0002-1677-6386<EMAIL_ADDRESS>url=https://www.bcglee.com/, ]

[orcid=0000-0002-9723-4262<EMAIL_ADDRESS>]

[orcid=0000-0001-7548-0932<EMAIL_ADDRESS>url=https://www.sarahhsalter.net/, ]

[orcid=0000-0001-5847-5242<EMAIL_ADDRESS>url=https://www.jim-
casey.com/, ]

# Navigating the Mise-en-Page: Interpretive Machine Learning Approaches to the
Visual Layouts of Multi-Ethnic Periodicals

Benjamin Charles Germain Lee∗ University of Washington  Joshua Ortiz Baco∗
Univeristy of Texas at Austin  Sarah H. Salter∗ Texas A&M University-Corpus
Christi  Jim Casey∗ Pennsylvania State University  $*$ = equal contribution

(2021)

###### Abstract

This paper presents a computational method of analysis that draws from machine
learning, library science, and literary studies to map the visual layouts of
multi-ethnic newspapers from the late 19th and early 20th century United
States. This work departs from prior approaches to newspapers that focus on
individual pieces of textual and visual content. Our method combines
Chronicling America’s MARC data and the Newspaper Navigator machine learning
dataset to identify the visual patterns of newspaper page layouts. By
analyzing high-dimensional visual similarity, we aim to better understand how
editors spoke and protested through the layout of their papers.

###### keywords:

machine learning periodicals newspapers editors document layout analysis
Chronicling America Newspaper Navigator MARC multi-ethnic press

##  1 Introduction

This paper presents a set of developing methods for analyzing Black, Latinx,
and other ethnic newspapers in the late nineteenth and early twentieth-century
United States. This historical period was marked by rising tides of racism and
xenophobia. Editing a newspaper dedicated to the progress of Black, Latinx, or
immigrant communities often entailed facing down forces of hate and violence.

Editors in the historical multi-ethnic press often chose to proceed
circumspectly. Rather than publish bold condemnations in headlines or
editorials, it might have been safer to rely on other affordances of
newspapers. It was possible to convey ideas using the paper’s very mise-en-
page, in between the lines (or columns) through editorial strategies such as
juxtaposing different content types, riffing on a given format, or
deliberately reprinting unusual or out-of-context material. In the multi-
ethnic press, editors often employed these strategies to render layout as its
own field of meaning-making: relations between items on the page could recast
or resignify seemingly distinct components of the newspaper. This feature of
multi-ethnic editorial communication makes it untenable to extract and analyze
isolated textual or visual items from these newspapers.

Computational methods make it possible to analyze the visual patterns of
multi-ethnic newspapers at scale. While prior research has approached these
questions through textual analysis, a pivot to computational visual analysis
compensates for the frequently noisy digital surrogates of surviving multi-
ethnic newspapers, too often poorly preserved.

Figure 1: A visualization of the front page of the May 22, 1917, issue of La
Prensa [1], annotated with bounding boxes of visual content derived from the
Newspaper Navigator dataset [2]. These bounding boxes are predictions made by
a machine learning model showing, in this case, the locations of headlines,
photographs, and maps (each bounding box includes the predicted class, as well
as a confidence score, in the top left corner). In this paper, we utilize
these bounding boxes to compare newspaper titles according to visual
similarity.

Our working method uses machine learning to map the layouts of hundreds of
thousands of newspaper pages, quantified by visual similarity, using the
Newspaper Navigator dataset (as pictured in Figure 1) [2]. Visualizing and
comparing layouts at scale provides a wayfinding function, helping to direct
archival researchers to distinctive and meaningful places in large collections
of digitized newspapers and cultural heritage collections writ large.

In this paper, we describe work in progress that brings machine learning,
library science, and literary studies into conversation. This work includes:
(1) a dataset of visual features from 16.3 million pages of historical
newspapers; (2) remediation in 207 catalog records for 309 ethnic editors; and
(3) modes for reading the practical languages of editing in specific
historical contexts. Prior research efforts blend together in our exploratory
machine learning prototypes for quantitatively and qualitatively measuring the
visual similarities of the formats of a newspaper page. These visual
similarities can help illuminate the ideas and activism that flowed through
the editorial craft and conduct of historical multi-ethnic newspapers.
Moreover, our work formulating a “similarity score” affirms the promise of
multi-disciplinary collaborations for improving access and analysis of larger-
scale digitized newspaper collections.

##  2 Related Prior Work on Developing Computational Approaches to Multi-
Ethnic Periodicals

Multiple fields are eager to explore the new possibilities for analyzing
digitized collections of historical newspapers. Complementary conversations in
computational, archival, and literary studies have begun to pursue those
possibilities in rich but largely disparate ways. Silos persist even as many
scholars make parallel use of the Library of Congress’s Chronicling America
and other national newspaper digital initiatives[3, 4, 5, 6, 7, 8, 9, 10, 11].

This paper builds on multi-ethnic press scholarship and archival research in
Chronicling America and other repositories. The prevailing view of political
activism in historical multi-ethnic newspapers has tended to focus on the
textual and visual content that drove demands for social justice [12, 13, 14,
15]. A complementary approach has focused on social histories, tracing the
networks of a paper’s contributors and readers [16, 17, 18, 19, 20]. What can
be easy to miss, however, is that editors can use the format of a paper itself
to voice critiques and protests. The formal nature of editors’ craft is
difficult to grasp in the meticulous approaches to archival research on a
case-by-case basis. Our work aims to develop experimental machine learning
applications to map the larger patterns of editorial practices and editing
conventions that span languages, communities, and eras.

Conversations in library and archival studies are exploring how best to
account for histories of race and ethnicity in metadata. Rather than being
reduced to a theme or topic within a collection, specialists have rethought
how race and ethnicity might work as organizing principles for collecting and
describing historical materials. Race and ethnicity as categories of identity
are historically contingent and socially constructed. That fluidity challenges
the rigors of knowledge organization systems that impose stable categories and
schemas. For some archivists, however, categories of race and identity
productively serve as a form of provenance to guide collection development and
access policies [21, 22, 23]. Many of these changes at the collection level
are being addressed in the metadata and, more recently, through linked open
data. Because the development of controlled vocabularies, such as the Library
of Congress Subject Headings, and Machine Readable Cataloging (MARC) largely
predate sensitivity to these histories in libraries and archives, metadata has
become a vital space for rediscovering the historical editors of multi-ethnic
newspapers [24].

At the same time, an overabundance of work with more “computationally robust,”
institutionally privileged white print collections has prompted preliminary
conversations regarding the directions of computational humanities research.
Indeed, from OCR engines that do not support indigenous languages or different
dialects, to entire subfields of natural language processing and textual
analysis that assume English as the default language of study,
algorithmically-mediated erasure continues to be uncovered and detailed by
scholars [25, 26, 27, 28, 29, 30, 31, 32, 33]. Many have advocated for the
development of computational approaches for those collections that remain
marginal. This paper concerns one such corrective for digitized periodicals:
utilizing machine learning to study visual content and page layout rather than
textual content, circumventing the use of OCR engines entirely.

Computational scholars have isolated visual content within periodicals in
order to produce datasets of isolated content. Humanistic scholars have
utilized these atomized items for interpretation. We are instead interested in
the visual layouts themselves and the relations they imply on the page. By
treating relations and patterns within layouts as forms of interpretive
content – reading them as a “visual grammar” of photographs, advertisements,
headlines, and beyond – we are in effect exploring a publication’s creative
history, local contexts, and circulation networks.

##  3 Methods

The work described in this paper began in part with the Newspaper Navigator
project, created as part of Library of Congress’s Innovator in Residence
program in 2019 with the goal of excavating visual content in Chronicling
America. Developed in partnership with LC Labs, the National Digital Newspaper
Program, and IT Design & Development at the Library of Congress, Newspaper
Navigator utilizes a machine learning model to extract seven visual features
from 16 million pages in Chronicling America: advertisements, photographs,
illustrations, maps, comics, editorial cartoons, and headlines [2]. This
machine learning model consists of an object detection model that had been
finetuned on thousands of bounding box annotations of visual content created
by volunteers as part of the Beyond Words crowdsourcing initiative launched by
LC Labs in 2017 [34]. The resulting extracted visual content was publicly
released as the Newspaper Navigator dataset, available in the public domain
for all to use [2].

Newspaper Navigator coincides with current trends in periodical studies. As
Sean Latham and Robert Scholes argue in an influential essay, “We have often
been too quick to see magazines [and other periodicals] as containers of
discrete bits of information rather than autonomous objects of study” [35]. In
practice, this means that periodical scholars tend to examine not only the
individual pieces of content in a given publication, but also its overall
construction, taking into account such architectural elements as a
publication’s formats, columns, and visual layout. The typical methods to
apply this approach require intensive reading of individual newspaper pages.
Researchers gradually develop an intuition about the general style and flavor
of a certain publication. A publication’s implicit patterns can be relatively
intangible because periodical formats are necessarily fixed—to preserve
continuity—and flexible—to accommodate new and changing scenarios. A paper
looks like itself almost always, but not always. Even beyond academic inquiry,
dedicated readers of any periodical often develop this kind of intuitive
mental map of a publication’s page designs and formats.

We saw the possibility of combining these two types of reading. On one hand,
Newspaper Navigator makes it possible to chart visual features across the
entirety of Chronicling America. On the other hand, the methods in periodical
studies encourage attention to the juxtaposition of visual features on a page.
Starting with these two methods, we developed a set of preliminary working
questions. How do the visual features on the page of a newspaper interlock and
mutually inform each other? If we look at all of those features holistically,
what kinds of styles and ideas can we detect over the lifespan of a single
newspaper? More broadly, what added perspective can we gain by using the
Newspaper Navigator to trace the sum of newspaper layouts over hundreds of
publications? Finally, can we establish a “baseline” for newspaper production
in the later 19th and early 20th century United States that would allow us to
see the critical departures made by multi-ethnic newspapers in the service of
social justice during moments of elevated racism or xenophobia?

###  3.1 Methods I: Adapting MARC 21 Metadata for Research

Adapting the Newspaper Navigator dataset to answer these questions required us
to collect additional data about the subsets of multi-ethnic newspapers within
Chronicling America. This information was partially extracted from the
bibliographic data of newspaper titles in the machine-readable cataloging
(MARC) format. All records for newspapers in Chronicling America follow the
cataloging standards of the Cooperative ONline SERials Program (CONSER) [36],
providing a somewhat normalized set of features for sorting newspapers.
Additionally, Chronicling America’s partnering institutions provide a “title
essay” that describes and contextualizes each newspaper. These unstructured
data has recently been the focus of efforts for metadata enhancement in
remediation projects geared to identifying ethnic editors [37].

The last step in collecting and organizing our data was combining MARC fields
and structured data generated from the text of title essays to find relevant
publications. Metadata from the title essays allowed us to select only ethnic
newspapers with identifiable editors. From this subset, we selected MARC
fields used by cataloguers to record title changes, race, ethnicity, and
languages associated with each publication. We initially hoped to find an
indicator of visual change within the metadata for newspapers. We considered
that serials with more than one language edition, language notes, or different
titles over their publication runs would feature changes to visual and layout
elements of newspapers to mark these relationships or changes.111The 780 and
785 fields indicate all preceding and succeeding titles for a publication,
respectively. The 546 field describes the language of a resource while the 765
indicates the original language of a translated publication and the 775
captures other language editions of a publication Title changes, for example,
are a cataloguing practice for linking preceding and succeeding titles of a
publication that has had a substantial name change, which we assumed would
also produce modifications to the graphic elements on the page. We were able
to parse these metadata fields where the change inherent in the seriality of
newspapers is recorded, which in turn provided parameters for capturing
difference, similarity, and linking in an otherwise static knowledge
organization system.

Figure 2: A heatmap of advertisements appearing on the front page of the
Opelousas Courier, 1870-1909. Darker regions on the heatmaps correspond to a
higher concentration of advertisement pages appearing in that region in
aggregate over a given year.

###  3.2 Methods 2: Quantifying and Visualizing the Mise-en-Page

The first visual output we experimented with were heat maps, created with
Newspaper Navigator data that encodes the positions of visual features on a
newspaper page. Inspired by PageOneX [38], these heat maps concisely summarize
visual content patterns in an interpretable fashion. Figure 2 shows the
locations and relationships of advertisements across four decades of The
Opelousas Courier,222The Opelousas Courier is described in Chronicling America
as a Louisiana multilingual newspaper, with content in French and English.
https://chroniclingamerica.loc.gov/lccn/sn83026389/marc/.  As an initial
experiment in visual layout representation, the heat maps enabled us to trace
the development of discrete visual features, otherwise impossible to quantify
by examining front pages by eye.

This approach facilitates several comparative opportunities. We rendered
multiple visual categories into heat maps to understand the interplay between
visual components in a single publication. We also utilized the heat map
graphics to compare different newspapers using the same visual categories. The
heat maps offered a useful preliminary system for studying newspaper layout as
historical information.

Heat maps added greater nuance to our historical questions about newspaper
layout. For example, what are the relationships between visual content,
editorial practices, and print technology? What are the roles of regional and
local context in predicting layout? How do we account for journalistic trends,
such as the rise of syndicated content or nationally-circulating
advertisements?

The heat maps showed the need for further exploration. In analyzing one
publication at a time, the heat maps reproduced the methods of archival
periodical studies, even in visually abstracted ways. While the MARC records
helped us to identify the titles we were most interested in rendering into
heat maps, this piecemeal approach to newspaper metadata and layouts proved
only a first step.

###  3.3 Methods 3: Toward a Visual Constellation of Multi-Ethnic Newspapers

Our next experiment attempted to scale up possible representations to refocus
on illustrating change and similarity across tens, hundreds, and even
thousands of publications. Our resulting exploratory model utilized the
quantified heat maps to define a metric over visual similarity. In particular,
the distance between two newspaper titles can be captured by the residuals
from subtracting two heatmaps from one another. As an initial approach, we
segmented each newspaper title’s run into individual years. For each title and
year, we then treated the seven resulting heat maps (one for advertisements,
one for headlines, etc.) collectively as a single, high-dimensional vector. In
order to interpret the clusters in this high-dimensional space, we utilized
T-SNE for dimensionality reduction [39].

Figure 3: A 2-dimensional map of visual similarity across all ethnic titles in
Chronicling America, filtered using [37]. Each point on the similarity map
represents the composite front page of a given newspaper title for a given
year, from 1890 to 1909 (individual years are labeled on the visualization).
The magnified cluster reveals that The Washington Bee, The Appeal, and The
Colored American are largely grouped together across years (individual front
pages from each title are shown in Figure 4). Notably, all three titles are
from the Black press. Figure 4: Sample front pages from The Appeal, The
Washington Bee, and The Colored American from 1891, 1896, and 1900,
respectively – all of which are title-year pairs appearing in the same cluster
in Figure 3 [40, 41, 42, 43, 44, 45]. All six front pages feature visual
content – namely, illustrations and photographs of individuals – prominently
in the center of the page.

This process resulted in exploratory constellations of newspaper titles, as
shown in Figure 3. In this figure, we also show a close-up of one such
“similarity cluster,” revealing that three Black newspapers (The Washington
Bee, The Appeal, and The Colored American) are largely grouped together across
years. In Figure 4, we show example front pages for these three titles from
years included within the cluster. All six pictured front pages feature visual
content on the center of the page. Notably, all six feature photographs and
illustrations of individuals, a trend that is even more apparent when browsing
these front pages across different issues. These images demonstrate that the
similarity clusters, generated using Newspaper Navigator’s rubric of seven
content classes, accurately reflect common visual patterns of layout in
distinct newspapers across the Black newspapers held in Chronicling America.
The shared qualities across each of the three newspapers pictured in Figure 4
are suggestive. They are roughly contemporaneous African American newspapers
that locate specific types of images near the center of the front page. These
results indicate the capacity of the method to produce a multifaceted research
question: what can we learn from the use of visual culture and portraiture in
the material pages of late nineteenth-century African American newspapers?

These constellations of similarity provide a means to uncover small networks
to focus on and interpret at a much more granular level. This is not a fully
quantitative representation, but rather an exploratory framework for
identifying and interpreting networks of related publications and the
editorial practices that shape them. Given our overarching interest in the
historical trajectories and meanings of newspaper formats, these
visualizations show great promise. The constellation map is a model for how to
convert raw bounding box information into interpretable provocations about
historical similarity. Similarity clusters introduce visual layouts as clues
to a shared context, inviting new research into the different historical
traditions, technological moments, or political causes that comprise the
histories of editorship.

##  4 Findings & Future Work

Our ongoing collaboration offers three central findings to date. First, this
work finds that full level CONSER MARC records provide useful data features in
machine-learning applications. Specifically, the most value comes from fields
that register change over time, such as preceding and succeeding titles, and
describe relationships between different editions of the same publication.
Further metadata enhancement focused on identifying ethnic editors will
significantly improve the interpretative value of visualizations and pattern
recognition from large-scale collections like Chronicling America.

Second, this collaboration reveals how computational notions of layout
similarity can be leveraged for humanistic inquiry. Using the machine
learning-constructed Newspaper Navigator dataset of visual content, we have
generated heat maps and constellation maps for exploring trends of visual
layout within and across newspaper titles, respectively. These visualizations
provide us with new visual grammars and affordances for excavating editorial
practices.

Third, this work required us to find a collaborative working process for
creating multidisciplinary methods of analysis that engage with machine
learning, library science, and the humanities. Rather than trying to blend or
transplant our respective research methods, we focused on creating a shared
vocabulary that could speak to current research agendas in each of the
relevant fields. This set of theoretical exchanges required iterative
conversations, frequently pedagogical in nature, to find a method that is not
merely extractive and static but generative and dynamic. This process will
help guide the future growth of conversations in computational periodicals
research.

Our roadmap for future work is guided by the provocations offered throughout
this paper. First, we plan to further utilize MARC data and the constellation
maps in concert with one another in order to inform our understanding of both.
The MARC data has the capacity to guide our navigation and understanding of
the constellation maps; conversely, the constellation maps provide mechanisms
for record enhancement by allowing scholars, catalogers, and librarians alike
to sift through incomplete legacy records of ethnic publications. We offer
this direction of future work in pursuit of uncovering editorship in ethnic
presses. We plan to refine our similarity metrics and visualizations from a
computational perspective in pursuit of this goal as well. The results of this
research will drive new inquiry into the longer and varied histories of
radical dissent and veiled protest in the multi-ethnic press.

###### Acknowledgements.

The authors would like to thank and acknowledge Dr. Molly Hardy, Senior
Program Officer at the National Endowment for the Humanities, along with
Georgia Higley and Vanessa Mitchell from the Serial and Government
Publications Division at the Library of Congress, for sharing their expertise.
This work would not be possible without Dr. Hardy’s input and advice.

## References

  * of Congress [1917] L. of Congress, La prensa, 22 may 1917. chronicling america: Historic american newspapers., 1917. URL: https://chroniclingamerica.loc.gov/lccn/sn83045395/1917-05-22/ed-1/seq-1/.
  * Lee et al. [2020] B. C. G. Lee, J. Mears, E. Jakeway, M. Ferriter, C. Adams, N. Yarasavage, D. Thomas, K. Zwaard, D. S. Weld, The Newspaper Navigator Dataset: Extracting Headlines and Visual Content from 16 Million Historic Newspaper Pages in Chronicling America, Association for Computing Machinery, New York, NY, USA, 2020, p. 3055–3062. URL: https://doi-org.offcampus.lib.washington.edu/10.1145/3340531.3412767.
  * of Congress & National Endowment for the Humanities [2021] L. of Congress & National Endowment for the Humanities, Chronicling America, 2021. URL: https://chroniclingamerica.loc.gov/.
  * Cassidy [2016] S. Cassidy,  Publishing the trove newspaper corpus,  in: Proceedings of the Tenth International Conference on Language Resources and Evaluation (LREC’16), European Language Resources Association (ELRA), Portorož, Slovenia, 2016, pp. 4520–4525. URL: https://www.aclweb.org/anthology/L16-1715.
  * Pekárek and Willems [2012] A. Pekárek, M. Willems,  The Europeana Newspapers – A Gateway to European Newspapers Online,  in: M. Ioannides, D. Fritsch, J. Leissner, R. Davies, F. Remondino, R. Caffo (Eds.), Progress in Cultural Heritage Preservation, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, pp. 654–659.
  * Willems and Atanassova [2015] M. Willems, R. Atanassova,  Europeana Newspapers: searching digitized historical newspapers from 23 European countries,  Insights 28 (2015) 51–56. URL: http://insights.uksg.org/articles/10.1629/uksg.218/. doi:10.1629/uksg.218, number: 1 Publisher: UKSG in association with Ubiquity Press.
  * Kreitz [2018] K. Kreitz,  American Alternatives: Participatory Futures of Print from New York City’s Nineteenth-Century Spanish-Language Press,  American Literary History 30 (2018) 677–702. doi:https://doi.org/10.1093/alh/ajy032.
  * del [NA] About delpher, N/A. URL: https://www.delpher.nl/nl/platform/pages/helpitems?title=wat+is+delpher.
  * bri [NA] About the british newspaper archive, N/A. URL: https://www.britishnewspaperarchive.co.uk/help/about.
  * our [NA] Ourdigitalworld: Digital newspapers, N/A. URL: https://ourdigitalworld.net/what-we-do/digital-newspapers/.
  * pap [NA] Papers past, https://natlib.govt.nz/collections/a-z/papers-past, N/A.
  * Nerone and Barnhurst [2003] J. Nerone, K. G. Barnhurst,  US newspaper types, the newsroom, and the division of labor, 1750–2000,  Journalism Studies 4 (2003) 435–449. URL: https://doi.org/10.1080/1461670032000136541. doi:10.1080/1461670032000136541, publisher: Routledge _eprint: https://doi.org/10.1080/1461670032000136541.
  * Cole [2020] J. L. Cole, How the Other Half Laughs: The Comic Sensibility in American Culture, 1895-1920, U Press of Mississippi, Jackson, 2020\.
  * Williams [2015] A. N. Williams,  Cultivating Black Visuality: The Controversy over Cartoons in the Indianapolis Freeman,  American Periodicals: A Journal of History & Criticism 25 (2015) 124–138. URL: http://muse.jhu.edu/article/601461. doi:10.1353/amp.2015.0032.
  * Stein [2007] S. A. Stein,  Making jews modern: The yiddish and ladino press in the russian and ottoman empires, the modern jewish experience,  International Journal of Middle East Studies 41 (2007). doi:10.1017/S0020743809090874.
  * Foster [2005] F. S. Foster,  A Narrative of the Interesting Origins and (Somewhat) Surprising Developments of African-American Print Culture,  American Literary History 17 (2005) 714–740. URL: http://muse.jhu.edu/journals/american_literary_history/v017/17.4foster.html, <p>Volume 17, Number 4, Winter 2005</p>.
  * Gardner [2015] E. Gardner, Black Print Unbound: The Christian Recorder, African American Literature, and Periodical Culture, Oxford University Press, 2015\. Google-Books-ID: 76n2CQAAQBAJ.
  * Rhodes [1998] J. Rhodes, Mary Ann Shadd Cary: the Black Press and Protest in the Nineteenth Century, Indiana University Press, Bloomington, IN, 1998.
  * Spires [2019] D. R. Spires, The practice of citizenship: Black politics and print culture in the early United States, University of Pennsylvania Press, Philadelphia, 2019.
  * Fagan [2016] B. Fagan, The Black Newspaper and the Chosen Nation, University of Georgia Press, Athens, GA, 2016.
  * Wurl [2005] J. Wurl,  Ethnicity as Provenance: in Search of Values and Principles for Documenting the Immigrant Experience,  Archival Issues 29 (2005) 65–76. URL: http://www.jstor.org/stable/41102095, publisher: Midwest Archives Conference.
  * Daniel [2014] D. Daniel,  Archival representations of immigration and ethnicity in North American history: from the ethnicization of archives to the archivization of ethnicity,  Archival Science 14 (2014) 169–203. URL: http://link.springer.com/10.1007/s10502-013-9209-6. doi:10.1007/s10502-013-9209-6.
  * Bastian [2014] J. A. Bastian,  Documenting Communities Through the Lens of Collective Memory,  in: D. Daniel, A. S. Levi (Eds.), Identity palimpsests: archiving ethnicity in the U.S. and Canada, number book 5 in Archives, archivists and society, Litwin Books, Sacramento, CA, 2014, pp. 15–34.
  * Ventura et al. [2019] G. B. Ventura, L. Gauthereau, C. Villarroel,  Recovering the US Hispanic Literary Heritage: A Case Study on US Latina/o Archives and Digital Humanities,  Preservation, Digital Technology & Culture (PDT&C) 48 (2019) 17–27. URL: https://www.degruyter.com/document/doi/10.1515/pdtc-2018-0031/html. doi:10.1515/pdtc-2018-0031.
  * Fagan [2016] B. Fagan,  Chronicling White America.,  American Periodicals 26 (2016) 10–13. URL: http://search.proquest.com/docview/2152559277/F95B44D158E64163PQ/13.
  * Alpert-Abrams [2016] H. Alpert-Abrams,  Machine Reading the Primeros Libros,  Digital Humanities Quarterly 010 (2016).
  * Joshi et al. [2020] P. Joshi, S. Santy, A. Budhiraja, K. Bali, M. Choudhury,  The state and fate of linguistic diversity and inclusion in the NLP world,  in: Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, Association for Computational Linguistics, Online, 2020, pp. 6282–6293. URL: https://www.aclweb.org/anthology/2020.acl-main.560. doi:10.18653/v1/2020.acl-main.560.
  * Yakel [2003] E. Yakel,  Archival representation,  Archival Science 3 (2003) 1–25. URL: https://doi.org/10.1007/BF02438926. doi:10.1007/BF02438926.
  * Hardy and DiCuirci [2019] M. Hardy, L. DiCuirci, Critical Cataloging and the Serials Archive: The Digital Making of “Mill Girls in Nineteenth-Century Print” \- Archive Journal, 2019. URL: http://www.archivejournal.net/?p=8073.
  * Williams [2019] L. Williams,  What computational archival science can learn from art history and material culture studies,  in: 2019 IEEE International Conference on Big Data (Big Data), 2019, pp. 3153–3155. doi:10.1109/BigData47090.2019.9006527.
  * Gustafson [2015] K. L. Gustafson,  Ethnic newspaper producers face archiving challenges,  Newspaper Research Journal 36 (2015) 314–327. URL: https://doi.org/10.1177/0739532915600744. doi:10.1177/0739532915600744.
  * Risam [2018] R. Risam, New digital worlds: postcolonial digital humanities in theory, praxis, and pedagogy, Northwestern University Press, Evanston, Illinois, 2018.
  * Noble [2018] S. U. Noble, Algorithms of Oppression: How Search Engines Reinforce Racism, NYU Press, 2018\. URL: http://www.jstor.org/stable/j.ctt1pwt9w5.
  * Labs [2017] L. Labs, Beyond words (“mark”), 2017\. URL: http://beyondwords.labs.loc.gov/#/mark.
  * Latham and Scholes [2006] S. Latham, R. Scholes,  The rise of periodical studies,  PMLA/Publications of the Modern Language Association of America 121 (2006) 517–531. doi:10.1632/003081206X129693.
  * con [NA] About CONSER - Program for Cooperative Cataloging (Library of Congress), N/A. URL: https://www.loc.gov/aba/pcc/conser/about/.
  * Ortiz Baco [2019] J. Ortiz Baco, Title Essays, Linked Data, and the Ethnic Press in Chronicling America, 2019\. URL: https://www.neh.gov/blog/title-essays-linked-data-and-ethnic-press-chronicling-america.
  * Costanza-Chock and Rey-Mazón [2016] S. Costanza-Chock, P. Rey-Mazón,  Pageonex: New approaches to newspaper front page analysis,  International Journal of Communication 10 (2016) 28.
  * van der Maaten and Hinton [2008] L. van der Maaten, G. Hinton,  Visualizing data using t-sne,  Journal of Machine Learning Research 9 (2008) 2579–2605. URL: http://jmlr.org/papers/v9/vandermaaten08a.html.
  * of Congress [1891a] L. of Congress, The appeal, 14 march 1891. chronicling america: Historic american newspapers., 1891a. URL: https://chroniclingamerica.loc.gov/lccn/sn83016810/1891-03-14/ed-1/seq-1/.
  * of Congress [1891b] L. of Congress, The appeal, 9 may 1891. chronicling america: Historic american newspapers., 1891b. URL: https://chroniclingamerica.loc.gov/lccn/sn83016810/1891-05-09/ed-1/seq-1/.
  * of Congress [1896a] L. of Congress, The washington bee, 29 feb. 1896\. chronicling america: Historic american newspapers., 1896a. URL: https://chroniclingamerica.loc.gov/lccn/sn84025891/1896-02-29/ed-1/seq-1/.
  * of Congress [1896b] L. of Congress, The washington bee, 29 aug. 1896\. chronicling america: Historic american newspapers., 1896b. URL: https://chroniclingamerica.loc.gov/lccn/sn84025891/1896-08-29/ed-1/seq-1/.
  * of Congress [1900a] L. of Congress, The colored american, 7 april 1900\. chronicling america: Historic american newspapers., 1900a. URL: https://chroniclingamerica.loc.gov/lccn/sn83027091/1900-04-07/ed-1/seq-1/.
  * of Congress [1900b] L. of Congress, The colored american, 8 sept. 1896\. chronicling america: Historic american newspapers., 1900b. URL: https://chroniclingamerica.loc.gov/lccn/sn83027091/1900-09-08/ed-1/seq-1/.


Nonparametric Hamiltonian Monte Carlo

Carol Makoxford
Fabian Zaiseroxford
Luke Ongoxford

oxfordDepartment of Computer Science, University of Oxford, United Kingdom

Carol<EMAIL_ADDRESS>
Inference algorithm, MCMC, Tree

Probabilistic programming uses programs to express generative models whose posterior probability is then computed by built-in inference engines. A challenging goal is to develop general purpose inference algorithms that work out-of-the-box for arbitrary programs in a universal probabilistic programming language (PPL). The densities defined by such programs, which may use stochastic branching and recursion, are (in general) nonparametric, in the sense that they correspond to models on an infinite-dimensional parameter space. However standard inference algorithms, such as the Hamiltonian Monte Carlo (HMC) algorithm, target distributions with a fixed number of parameters. This paper introduces the Nonparametric Hamiltonian Monte Carlo (NP-HMC) algorithm which generalises HMC to nonparametric models. Inputs to NP-HMC are a new class of measurable functions called “tree representable”, which serve as a language-independent representation of the density functions of probabilistic programs in a universal PPL. We provide a correctness proof of NP-HMC, and empirically demonstrate significant performance improvements over existing approaches on several nonparametric examples.

1. Use $\dif$ (Roman) in integrals

2. Changed “bold blue path” to “blue path”

3. Nonparametric (not non-parametric)

4. initialise / formalise (not initialize / formalize) [British English, for consistency]

5. Use  instead of  if possible, for consistency

6. Use  for referencing Assumption #1.

§ INTRODUCTION

Probabilistic programming languages (PPLs)

are tools to study real-world random phenomena with observed data, by specifying probabilistic models using a mix of deterministic and stochastic constructs.

Unlike other representations of probabilistic models, PPLs benefit
from the simplicity and expressiveness of programming languages and makes probabilistic modelling clearer and more widely applicable.

PPLs usually come with built-in algorithms to perform Bayesian inference on such models automatically.
This frees

from the mathematically involved task of designing custom inference algorithms for their models.

[lo]Probabilistic programming is a general purpose means of expressing probabilistic models as programs, and automatically performing Bayesian inference.
Probabilistic programming systems enable data scientists and domain experts to focus on designing good models;
the task of developing efficient inference engines can be left to experts in Bayesian statistics, machine learning and programming languages.
To realise the full potential of probabilistic programming, it is essential to automate the inference of latent variables in the model, conditioned on the observed data.

[lo]Church <cit.> introduced universal probabilistic programming, the idea of writing probabilistic models in a Turing-complete functional programming language.
[lo]Typically containing only a handful of basic programming constructs such as branching and recursion, universal probabilistic programming languages (PPLs) can nonetheless specify all computable probabilistic models <cit.>.

In particular, nonparametric models—models with an unbounded number of random variables—can be described naturally in universal PPLs using recursion.
These include probabilistic models with an unknown number of components, like
Bayesian nonparametric models <cit.>,
variable selection in regression <cit.>,
signal processing <cit.>;
and models that are defined on infinite-dimensional spaces, such as probabilistic context free grammars <cit.>,

birth-death models of evolution <cit.> and
statistical phylogenetics <cit.>.

Examples of practical universal PPL include Anglican <cit.>, Venture <cit.>, Web PPL <cit.>, Hakaru <cit.>, Pyro <cit.>, Turing <cit.> and Gen <cit.>.

[fz]While universal PPLs offer a lot of flexibility when specifying models, this makes
designing efficient general-purpose Markov chain Monte Carlo (MCMC)
Why MCMC? Isn't general-purpose enough? inference algorithms for
[cm]The flexibility offered by universal PPL makes

designing efficient general-purpose Markov chain Monte Carlo (MCMC) Why MCMC? Isn't general-purpose enough? Because we are not comparing with DCC
inference algorithms for

[fz]them very challenging.
However, because universal PPLs are expressively complete, it is a challenging problem to design and implement general purpose inference engines for them.
The parameter space of a nonparametric model is a union of spaces of varying dimensions.

[fz]To approximate the posterior via an Markov chain Monte Carlo (MCMC) algorithm, the transition kernel will have to efficiently switch between a potentially unbounded number of configurations of different dimensions.

This difficulty explains why there are so few general purpose MCMC algorithms for universal PPLs <cit.>.
[lo]We believe it is also the reason why these algorithms struggle with nonparametric models, as we show in <ref>.
A case in point is the widely used universal PPL Pyro.
Even though it allows the specification of nonparametric models, its HMC and No-U-Turn Sampler <cit.> inference engines do not support them reliably:

[fz]in one of our benchmark tests, they produced a wrong posterior (<ref>).

In this paper, we introduce the Nonparametric Hamiltonian Monte Carlo (NP-HMC) algorithm, which generalises the Hamiltonian Monte Carlo (HMC) algorithm <cit.> to nonparametric models.
The input to NP-HMC is what we call a tree representable (TR) function, which is a large class of measurable functions of type [fz]$\btra: \bigcup_{n \in \Nat} \Real^n \to \Real_{\ge 0}$,
designed to be a language-independent representation for the density functions of programs written in any universal PPL.
The parameter space of the standard HMC algorithm is $\Real^n$, a Euclidean space of a fixed dimension.
By contrast, the parameter space of NP-HMC is $\bigcup_{n \in \Nat}\Real^n$.
[lo]The key innovation of NP-HMC is a method by which the dimension of the configuration of the current sample is incremented lazily, while preserving the efficacy of HMC by keeping the Hamiltonian approximately invariant.
We prove that NP-HMC is correct, i.e., the induced Markov chain converges to the posterior distribution.

To evaluate the practical utility of NP-HMC, we compare an implementation of the algorithm against existing out-of-the-box MCMC inference algorithms on several challenging models with an unbounded number of random variables.
Our results suggest that NP-HMC is applicable to a large class of probabilistic programs written in universal PPLs,
offering significantly better performance than existing algorithms.

Motivated by the effectiveness of the Hamiltonian Monte Carlo (HMC) inference algorithm <cit.>,
we address this shortcoming by proposing the Nonparametric Hamiltonian Monte Carlo (NP-HMC) algorithm and show that it is a (provably) correct and practically effective inference algorithm for universal PPLs.
To set the scene, we introduce the notion of tree representable (TR) functions

[fz]as a language-[cm]independent

abstraction for the density functions of programs written in universal PPLs.
We then demonstrate how NP-HMC generalises the standard HMC algorithm to nonparametric models specified by TR functions:

[fz]it lazily extends or should we say dynamically instead of lazily? the dimension of the parameter spacewhat does "according to the current sample" mean?, all while preserving the effectiveness of HMC.
[fz]We prove that NP-HMC is correct, i.e. that its Markov chain converges to the posterior distribution.
Furthermore, we show that [fz]our adoption of HMC is conceptually generic by presenting the nonparametric extensions of [fz]two HMC variants: reflective/refractive HMC <cit.> and discontinuous HMC <cit.>.
To evaluate the practical utility of NP-HMC, we compare an implementation of it against existing

[fz]out-of-the-boxnot just MCMC inference algorithms on challenging models with an unbounded number of random variables.
Our results suggest that NP-HMC can be effectively applied to
[fz]a large class of probabilistic programs written in universal PPLs,
offering better performance than existing algorithms.
Claiming that it is applicable to arbitrary universal probabilistic programs goes a bit far, IMO, considering how much we struggled to find working examples at the start.
Built-in inference algorithms are powerful tools in probabilistic programming language (PPL) frameworks, such as Stan <cit.> and Pyro <cit.>.

Such inference algorithms include Lightweight Metropolis-Hastings <cit.>, Reversible Jump MCMC <cit.>,
but one of the most popular is the Hamiltonian Monte Carlo (HMC) algorithm <cit.>.
It leverages the recent advances of automated gradient computation to efficiently sample from target distributions with densities of type $\Real^n \to \pReal$.
However, HMC only works for a fixed dimension of the parameter space; it cannot transition between traces of different [lo]and unbounded lengths.
In this paper, we address this shortcoming.

LAFI reviewer 1: What is the main technical shortcoming of HMC that is solved by NP-HMC?

We should say that this work is motivated by probabilistic programming.
Ultimately, we want to develop a general inference algorithm that works for all programs of a universal PPL.

To apply HMC on universal probabilistic programs,
we first introduce the notion of tree representable (TR) function which identifies the class of measurable functions of type $\bigcup_{n\in\Nat}\Real^n \to \pReal$ that arise as the density function (in the sense of trace semantics <cit.>) of programs written in universal PPLs, allowing recursion and (stochastic) branching, such as <cit.>.
Then, we present the Nonparametric Hamiltonian Monte Carlo (NP-HMC) algorithm,
which, given a , uses Hamiltonian dynamics to generate a Markov chain that converges to the distribution induced by it.
We present a correctness proof of NP-HMC and some preliminary experimental results based on a prototype implementation.

[To be continued.]

We write $\seqa$ to mean a (possibly infinite) real-valued sequence;
$\seqrange{1}{i}$ [lo]the prefix of $\vec \seqa$ consisting of the first $i$ coordinates;
$\seqi{i}$ the $i$-th coordinate of $\seqa$; and
$\len{\seqa}$ the length of $\seqa$.
We write sequence as lists,
such as $[3.6,1.0,3,55, -4.2]$, and

the concatenation of sequences $\seqa$ and $\seqa'$ as $\seqa \concat \seqa'$.

We write $\Borel_n$ for the Borel $\sigma$-algebra of $\Real^n$;

$\Gau_n$ for the standard $n$-dimensional normal distribution with mean $\boldsymbol 0$ and covariance $\idmat$;
[cm]$\pdfGau_n(\vec{x} \mid \vec{\mu}, \vec{\Sigma})$

for the density of $\vec{x} \in \Real^n$ in the $n$-dimensional normal distribution with mean $\vec{\mu}$ and covariance $\vec{\Sigma}$.
For brevity we write $\pdfGau_n(\vec{x})$ for $\pdfGau_n(\vec{x} \mid \vec{0}, \vec{\idmat})$ and $\pdfGau$ for $\pdfGau_1$.

Can we inline some of the definitions in this paragraph to where they're used? It seems like some are only used once.

For any $\pReal$-valued function
[lo]$f: \domain{f} \to \pReal$,
we write $\support{f} := \inv{f}(\Real_{>0})$ for the support of $f$; and

$\nsupport{f}{n} := \support{f} \cap \Real^n$ for the support of $f$ in $\Real^n$.
We say $x \in X$ is $f$-supported if $x \in \support{f}$.

§ TREE REPRESENTABLE FUNCTIONS

Conventional HMC [fz]samples from a distribution with a density function $\btra: \Real^n \to \Real_{\ge 0}$ where the dimension of the target (parameter) space $\Real^n$ is fixed.
However this is too restrictive for probabilistic programs, because—with branching and recursion—the target space has a variable, even unbounded, number of dimensions.

[caption=A simple probabilistic program.,label=intro-program,captionpos=b,float]
q = sample(normal(0, 1))
sum = 0
while sum < q:
sum += sample(normal(0, 1))
observe(sum, normal(q, 1))

Consider the probabilistic program in <ref> where sample(normal(0, 1)) denotes sampling from the standard normal distribution.

The dimension of the target space, i.e. the number of samples drawn, can vary from run to run because of the branching behaviour of the while loop.
(Recall that by trace, we mean the sequence of samples drawn in the course of a particular run, one for each random primitive encountered.)
[fz]The density function then records the weight of this trace, computed by multiplying the probability densities of all sampled values and the likelihoods of all observations during the run.
For the above program, we could have a trace $[0.3, 0.5] \in \Real^2$ of length 2 or a trace $[1.0, 0.5, 0.5] \in \Real^3$ of length 3, and so on.
Hence we have to consider density functions of type $\btra: \bigcup_{n\in \Nat} \Real^n \to \Real_{\ge 0}$.
However, not every such function makes sense as the density of a probabilistic program.
For example, if $\btra([1.0, 0.5, 0.5]) > 0$ then the program can execute successfully with the trace $[1.0, 0.5, 0.5]$, but not with $[1.0, 0.5]$ or any other prefix.
[fz]In other words, no proper prefix of $[1.0, 0.5, 0.5]$ is in $\support{\btra}$.

Thus we set our target space to be the measure space of traces
$\traces := \bigcup_{n\in\Nat} \Real^n$ equipped with the standard disjoint union $\sigma$-algebra
$\Sigma_\traces := \set{\bigcup_{n\in\Nat} U_n \mid U_n \in \Borel_n }$, with measure given by summing the respective (higher-dimensional) normals $\tmeasure(\bigcup_{n\in\Nat} U_n) := \sum_{n\in\Nat} \Gau_n(U_n)$.

We consider density functions that are measurable functions $\btra: \traces \to \pReal$ satisfying the prefix property:

whenever $\vec{q} \in \nsupport{\btra}{n}$ then for all $k < n$,
we have $\vec{q}^{1\dots k} \not\in \nsupport{\btra}{k}$.

We call them tree representable (TR) functions because any such function $\btra$ can be represented as a possibly infinite but finitely branching tree, which we call program tree.

This is exemplified in <ref> (left),
a circular node denotes an element of the input of type $\Real$;
a rectangular node gives the condition for
$\vec{q} \in \nsupport{\btra}{n}$
(with the left, but not the right, child satisfying the condition); and
a leaf node gives the result of the function on that branch.
Any branch (i.e. path from root to leaf) in a program tree of $w$ represents a set of finite sequences $[q_1,\dots,q_n]$ in $\support{\btra}$.
In fact, every program tree of a TR function $w$ specifies a countable partition of $\support{\btra}$ [fz]via its branches.
The prefix property guarantees that [fz]for each TR function $w$, there are program trees of the form in <ref> (left) representing $w$.

- Why call it “binary tree” representable (there seems no connection with the defining axiom.)? Explain the name  by presenting the “binary tree program” that computes a  $\btra$.
BTW we will first need to define $w_{\le n}$ before presenting the program.

- Is  a naturally large class of functions? Yes: densities denoted by SPCF programs (with standard normal distributions) are tree representable. (Give / sketch the algorithm that compiles SPCF programs into .)
From a PL perspective,  may be viewed as a PL-neutral IR (or intermediate representation) for our new inference algorithm, NP-HMC.

- Why do we need / how do we use  in NP-HMC? Our approach (the algorithms) assumes that the input density function is organised into a family of potential energies. Any  can be so organised.

path/.style=draw=blue, line width=3pt
not path/.style=draw=black, line width=0.5pt
input/.style=circle,draw=black, line width=0.5pt,minimum size=12pt,font=,inner sep=2pt
condition/.style=rectangle,draw=black, line width=0.5pt, minimum height=0.1pt,font=,inner sep=2pt

[scale=.55,sibling distance = 5em]
[input] (q1) $q_1$
node [condition] (c1) $[q_1] \overset{?}{\in} \nsupport{\btra}{1}$
child node [leaf] (r1) $\btra([q_1])$
node [input] (q2) $q_2$
node [condition] (c2) $[q_1,q_2] \overset{?}{\in} \nsupport{\btra}{2}$
child node [leaf,xshift=-1em] (r2) $\btra([q_1,q_2])$
node [input] (q3) $q_3$
node [condition] (c3)
$[q_1,q_2,q_3] \overset{?}{\in} \nsupport{\btra}{3}$
child node [leaf,xshift=-1em] (r3) $\btra([q_1,q_2,q_3])$
child node $\vdots$

at ($(c1) + (220:1.3cm) $) [decision] yes;
at ($(c1) + (315:1.1cm) $) [decision] no;
at ($(c2) + (215:1.4cm) $) [decision] yes;
at ($(c2) + (315:1.1cm) $) [decision] no;
at ($(c3) + (215:1.4cm) $) [decision] yes;
at ($(c3) + (315:1.1cm) $) [decision] no;

[scale=.55,sibling distance = 5em]
[input] (q1) $q_1$
node [condition] (c1) $q_1 \leq 0$
edge from parent [path]
child node [leaf,xshift=-1em] (r1) $\pdfGau(0 \mid q_1,1)$
edge from parent [not path]
node [input] (q2) $q_2$
node [condition] (c2) $0<q_1\leq q_2$
child node [leaf,xshift=-1.25em] (r2) $\pdfGau(q_2 \mid q_1,1)$
node [input] (q3) $q_3$
edge from parent [not path]
node [condition, text width=2.1cm] (c3)
$0 < q_1$ and  $q_2 < q_1 \leq q_2+q_3$
child node [leaf,xshift=-2em] (r3) $\pdfGau(q_2+q_3 \mid q_1,1)$
child node $\vdots$

[position, right of=q1, node distance=17pt] $0.3$;
[position, right of=q2, node distance=17pt] $0.5$;

(left) Generic  $\btra$; (right)  $\btra$ for the probabilistic program in <ref>.

We target s as densities for our new sampler NP-HMC because they are a naturally large class of functions.
In particular, every program of a universal PPL has a density function that is tree representable.[ With (additional) suitable assumptions about the computability of $w$, we can view any such tree as the abstract syntax tree of a program that computes $w$, but with any recursion unravelled (so that the tree is potentially infinite).]

<ref> in <ref>
the SM
for a formal account.)

For instance, the program in <ref> has density[Notice that, even though the program samples from a normal distribution, $\btra$ does not factor in Gaussian densities [fz]from those sample statements, just the observe statement, since they are already accounted for by $\measure{\traces}$.]
$\btra$ (.r.t. the [cm]stock measure $\tmeasure$) given by

\[
{\btra}{(\vec{q})} :=
\begin{cases}
{\pdfGau(\sum_{i=2}^{n} \vec{q}_i \mid \vec{q}_1,1)}
\begin{matrix}
\text{if }
\forall k < \len{\vec{q}} \text{\hspace{5em}} \\[-0.1em]
{\sum_{i=2}^{k} \vec{q}_i < \vec{q}_1 \leq \sum_{i=2}^{n} \vec{q}_i},
\end{matrix}
\\
0 & \text{otherwise}
\end{cases}
\]

which is TR and [lo]it has a program tree as depicted in <ref> (right).
Notice that every element in the support of $\btra$ belongs to a branch in this tree: for example, the trace $[0.3,0.5]$ belongs to the blue branch in <ref> (right).

As we will explain in <ref>, the prefix property (satisfied by TR functions) is

essential for the correctness of NP-HMC.

§ NONPARAMETRIC HMC

Nonparametric Hamiltonian Monte Carlo (NP-HMC) (<ref>) is a MCMC algorithm that, given a The qualification “(computable)” may prompt questions. In view of footnote <ref>, I think we can / should drop it.  $\btra$, iteratively proposes [lo]a new sample $\seqa \in \bigcup_{n \in \Nat}\Real^n$ and accepts it with a suitable Hastings acceptance probability,

[lo]such that the invariant distribution of the induced Markov chain is

\[ \tdist:{A}\mapsto \frac{1}{Z} \shortint{A}{\btra}{\tmeasure}
\]

with normalising constant $Z := \shortint{\traces}{\btra}{\tmeasure}$.
As the name suggests, NP-HMC is a generalisation of the standard HMC algorithm (<ref>.i) to nonparametric models, in the form of TR functions whose support is a subspace of $\traces$ of unbounded dimension.

In this section, we first explain our generalised algorithm, using a version (<ref>) that is geared towards conceptual clarity, and defer discussions of more efficient variants.

axis lines = left,
xlabel = $q_1$,
ylabel = Energy,
xmin=-4, xmax=4,
every axis plot/.append style=very thick,
grid style=line width=.1pt, draw=gray!10,
major grid style=line width=.2pt,draw=gray!50,

$-\log(\pdfGau(0 \mid q_1,1))\cdot[q_1 \leq 0]$

[label=45:[-3.1],circle,fill,inner sep=2pt] at (axis cs:-3.1,5.72393853) ;
[label=45:[-2.37],circle,fill,inner sep=2pt] at (axis cs:-2.37,3.72738853) ;
[label=45:[-0.86],circle,fill,inner sep=2pt] at (axis cs:-0.86,1.28873853) ;

Time    0    1    2    3

$\vec{q}$    [-3.1]    [-2.37]    [-0.86]    [1.15]

$\vec{p}$    [1.2]    [3.29]    [3.94]    [5.26]

The Hamiltonian dynamics of a particle on the surface $-\log(\pdfGau(0 \mid q_1,1))\cdot[q_1 \leq 0]$ on $\Real$.

axis lines = left,
xlabel = $q_1$,
ylabel = $q_2$,
zlabel = Energy,
xmin=-4, xmax=4,
ymin=-4, ymax=4,
legend style=at=(1.05,1.05),
every axis plot/.append style=thick,
grid style=line width=.1pt, draw=gray!10,
major grid style=line width=.2pt,draw=gray!50,
legend cell align=left,

3 [
y domain=-4:4,
$-\log(\pdfGau(0 \mid q_1,1))\cdot[q_1 \leq 0]$

3 [
y domain=0:4,
restrict z to domain=0:8,
](y > x) ? -ln(1/sqrt(2*pi)*exp(-(1/2)*(y-x)^2)):-1;
$-\log(\pdfGau(q_2 \mid q_1,1)) \cdot[0 < q_1 \leq q_2] $

[label=[fill=white, fill opacity=0.75, text opacity=1]90:[-3.1, -1.61],circle,fill,inner sep=2pt] at (axis cs:-3.1,-1.61,5.72393853) ;
[label=[fill=white, fill opacity=0.75, text opacity=1]90:[-2.37, -0.39],circle,fill,inner sep=2pt] at (axis cs:-2.37,-0.39,3.72738853) ;
[label=[fill=white, fill opacity=0.75, text opacity=1]90:[-0.86, 0.82],circle,fill,inner sep=2pt] at (axis cs:-0.86,0.82,1.28873853) ;

[label=[fill=white, fill opacity=0.75, text opacity=1]45:[1.15, 2.04],circle,fill,inner sep=2pt] at (axis cs:1.15,2.04,1.31321053) ;

Time    0    1    2    3

$\vec{q}$    [-3.1, -1.61]   [-2.37, -0.39]   [-0.86, 0.82]   [1.15, 2.04]

$\vec{p}$    [1.2, 3.04]    [3.29, 3.04]    [3.94, 3.04]    [5.26, 3.04]

The Hamiltonian dynamics of a particle on the updated surface $-\log(\pdfGau(0 \mid q_1,1))\cdot[q_1 \leq 0] -\log(\pdfGau(q_2 \mid q_1,1)) \cdot[0 < q_1 \leq q_2] $ on $\Real^2$.

[fz]We assume basic familiarity with the HMC algorithm; see <cit.> for the intuition behind it and <cit.> for details.
Like the HMC algorithm, NP-HMC views a sample $\seqa$ as a particle at position $\seqa$, with a randomly initialised momentum $\vec p$, moving on a frictionless surface derived from the density function $\btra$.
The key innovation lies in our treatment when the particle moves beyond the surface, i.e. outside the support of the density function.
This procedure is called $\extend$ (<ref>), which we will now illustrate.

Let's trace the movement of a particle
at position $\vec{q}=[-3.1]$ with a randomly chosen momentum $\vec{p} = [1.2]$,
on the [lo]surface (a line in 1D) determined by the  $\btra$ in <ref> (right).
The first two steps taken by the particle are simulated according to the Hamiltonian dynamics on the surface $-\log(\pdfGau(0 \mid q_1,1))\cdot[q_1 \leq 0]$, which is derived from the restriction of $\btra$ to $\Real$.
The positions on the surface and states[As in HMC, a state of the NP-HMC algorithm is a position-momentum pair $(\vec{q}, \vec{p})$ with $|\vec{q}|=|\vec{p}|$; but unlike HMC, $\vec{q}, \vec{p} \in \traces$.] of the particle at each step are given in <ref>.

At the third time step, the particle is at the position $[1.15]$, which is no longer on the surface.
To search for a suitable state, NP-HMC increments the dimension of the [lo]current surface (line) as follows.

First, the surface is extended to
$-\log(\pdfGau(0 \mid q_1,1))\cdot[q_1 \leq 0] -\log(\pdfGau(q_2 \mid q_1,1)) \cdot[0 < q_1 \leq q_2]$,
which is derived from the sum of the respective restrictions of $\btra$ to $\Real$ and to $\Real^2$,
as depicted in <ref>.
Since $\btra$ satisfies the prefix property, the respective supports of these restrictions,
namely $\set{[q_1,q_2]\in\Real^2 \mid q_1 \leq 0}$ and $\set{[q_1,q_2]\in\Real^2 \mid 0 < q_1 \leq q_2}$, are disjoint; and
hence the states of the particle on the previous surface $-\log(\pdfGau(0 \mid q_1,1))\cdot[q_1 \leq 0]$ can be reused on the updated surface.
[lo]Notice that the respective first coordinates of the particle's positions on the surface in <ref> are identical to the particle's positions on the surface (line) in <ref>.

Next, [lo]the initial state $(\vec{q} =[-3.1],\vec{p}=[1.2])$ is extended by appending a randomly chosen value to both the position and momentum components, so that the particle is positioned on the updated surface with an initial momentum.
In our example, $-1.61$ and $3.04$ are sampled and the initial state of the particle becomes $(\vec{q} =[-3.1,-1.61],\vec{p}=[1.2,3.04])$ which is located on the surface as shown in <ref>.
The states at times $1$, $2$ and $3$ are updated accordingly and are given in the table in <ref>.
Notice that the particle at time $3$ is now positioned on the updated surface, and hence Hamiltonian dynamics can resume.
The rest of this section is devoted to [lo]formalising our algorithm.

Henceforth we assume that the input  $\btra$ satisfies the following:

Integrability, i.e., the normalising constant $Z < \infty$.
(Otherwise, the inference problem would be ill-defined.)

The function $\btra$ is almost everywhere continuously differentiable on $\traces$.
(Since Hamiltonian dynamics exploits the derivative of $\btra$ in simulating the position of a particle,
this ensures that the derivative exists “often enough” for the Hamiltonian integrator to be used.)

For almost every infinite real-valued sequence $\seqa$,
there is a $k$ such that $\btra$ is positive on $\seqrange{1}{k}$.
(This ensures the $\extend$ subroutine in <ref> terminates almost surely.)

* NP-HMC is a generalisation of HMC to nonparametric models. Precisely, NP-HMC specialises to standard HMC (with the leapfrog integrator) in case the [fz]density function $w$ satisfies [cm]$\domain{w} = \Real^n$ for some $n$,

in addition to Assumptions ass:11, ass:22 and ass:33.

* All closed integrable almost surely terminating programs of a universal PPL[[cm]An almost surely terminating program (as defined in
<ref>the SM)
almost never observes a value with zero probability density.]
induce densities that satisfy Assumptions ass:11, ass:22 and ass:33.
[lo](See <ref> and <ref> the Supplementary Materials (SM) for an account.)

I don't think (3) holds for arbitrary SPCF programs. E.g. consider the program . I think it only holds if we additionally almost never score with 0.
Then (3) follows directly from almost sure termination.
Thanks for this. You're right and it does not hold for all SPCF programs. However, I think what you suggested is not enough, as we can still score a constant function at 0. It might be better/easier to change the small-step reduction of scoring to the following.
$\config{\Score{\PCF{r}}}{w}{s} \red \MyCase{\config{\PCF{r}}{w\cdot r}{s}}{r > 0}{\mathsf{fail}}$
In which case, almost surely terminating should be enough. Do you agree?

I agree and I think we mean the same thing. I meant that
I don't think we need to change the reduction rule for this though. An informal sentence like the emphasized one above should be enough.

The surfaces on which the particle is positioned are derived from [lo]a sum of appropriate restrictions of the input  $\btra$, defined as follows.
The $n$-th truncation $\trunc{n}:\Real^n \to \pReal$ is
$\seqa \mapsto \sum_{k=1}^n \btra(\seqrange{1}{k})$,
which returns the cumulative sum of the weight on the prefixes of an $n$-dimensional trace $\seqa$.
Thanks to the prefix property, for each $\vec q$, at most one summand is non-zero.
So any real-valued $\trunc{n}$-supported sequence $\seqa\in\Real^n$ has a prefix in the support of $\btra$;
and any $\btra$-supported sequence of length $n$ is also $\trunc{n}$-supported.

We define a family $U = \set{U_n}_{n\in\Nat}$ of potential energies
where each $U_n:\Real^n \partialto \pReal$ is a partial function defined as $U_n := -\log \trunc{n}$
with domain $\domain{U_n} := \support{\trunc{n}}$.
[cm]These are the surfaces on which the particle is positioned.

Proposal step
The nonparametric (NP) integrator $\BTint$ (<ref>) proposes a state by simulating the Hamiltonian motion of a particle at position $\seqa \in \Real^n$, with potential energy $U_n(\seqa)$ and a randomly chosen momentum $\vec{p} \in \Real^n$.[The Hamiltonian motion is almost always defined, by Ass. ass:22.]
The simulation runs $L$ discrete update steps [fz](also called leapfrog steps) of size $\epsilon > 0$, or until
the particle leaves the domain of $U_n$ (i.e. the support of $\trunc{n}$).
At that moment, the simulation stops and
NP-HMC “extends” the state $(\seqa,\vec{p})$ via the $\extend$ subroutine (<ref>), until the position of the particle falls into the domain of some potential energy [lo](i.e. support of some higher dimensional truncation).
[This happens almost surely, thanks to Ass. ass:33.]
Once the [fz]extended position of the particle is settled, simulation resumes.
[fz]If no extension is necessary, the behaviour is the same as that of the standard HMC leapfrog integrator.

The heart of NP-HMC is the $\extend$ subroutine (<ref>).

$\extend$ is called after
$i$ position steps and $i-1/2$ momentum steps are completed, i.e., at time $t = i \, \epsilon$.
If $\seqa$ is in the domain of $U_{\len{\seqa}}$, $\extend$ leaves the state unchanged;
otherwise [fz]$\seqa$ is not long enough and the while loop extends it as follows.

* Sample a pair $(x_0,y_0)$ of real numbers from the standard normal distribution respectively.

* Trace the motion of a particle with constant potential energy $1$
for $i$ position and [lo]$i-1/2$ momentum steps,
starting from $(x_0,y_0)$,
to obtain $(x,y)$.
Notice that the momentum update is simply the identity, hence we only need to consider $i$ position updates which takes $x$ to $x_0 + t \, y_0$.

* Append
$(x_0,y_0)$ to the initial state $(\vec{q_0},\vec{p_0})$ and
$(x,y)$ to the current state $(\seqa,\vec{p})$.

Thus the length of the position $\seqa$ is incremented,
and by [lo]Assumption ass:33, this loop terminates almost surely at a position $\seqa$ in the domain of the potential energy of dimension $|\seqa|$.

Putting them together
A single NP-HMC iteration, as shown in <ref>, produces a proposed sample $\vec{q}$ from the current sample $\vec{q_0}$, by applying the NP integrator $\BTint$ to the state $(\vec{q_0},\vec{p_0})$ with a freshly sampled momentum $\vec{p_0}$.
The proposed state $(\vec{q},\vec{p})$ is then accepted with probability given by the Hastings acceptance ratio.

The prefix property of the input  $\btra$ plays an important role in the NP-HMC inference algorithm.
If <ref> returns $\vec q$ as the next sample,
then for any extension $\vec q'$ of $\vec q$, we have $w(\vec q') = 0$, and so, all such $q'$ are irrelevant for inference ([fz]because $U_{|\vec q'|}(\vec{q'}) = U_{|\vec q|}(\vec{q})$).
If the prefix property weren't satisfied,
the algorithm would fail to account for the weight on such $\vec q'$.

Is this remark really helpful or can it be deleted?

Other extensions and efficiency considerations

We have presented a version of NP-HMC that eschews runtime efficiency in favour of clarity.

An advantage of such a presentation
[lo](in deliberately purified form)

is that it becomes easy to see that the same method is just as applicable to such HMC variants as
reflective/refractive HMC <cit.> and
discontinuous HMC <cit.>; we call the respective extensions NP-RHMC and NP-DHMC.
For details, see
the SM

I changed this.
Before we said we implemented NP-RHMC, which is only partially true.
I implemented it, but only for very simple branches of the form .
I barely tested this implementation.
The original RHMC paper implements it for general affine conditions, but this is difficult to implement and still insufficient for our experiments anyway (cf. stick-breaking).
Since NP-DHMC does not require computing intersections at all, I think it is just the superior algorithm.
Nevertheless, we should not claim to have fully implemented NP-RHMC.

Several efficiency improvements to NP-HMC are possible.
If the density function is given by a probabilistic program, one can interleave its execution with $\extend$, by gradually extending $\vec q$ at every encountered sample statement.
Similarly, the truncations $w_{\le n}$ don't require an expensive summation to compute.
For details, see
the SM

Our implementation (<ref>) also improves the $\extend$ function (<ref>): it not only extends a trace $\vec q$ if necessary, it also trims it to the unique prefix $\vec q'$ of $\vec q$ with positive $w(\vec q')$.
This version works better in our experiments.
For details, see
the SM

Because of space restriction, we may need to move this to the implementation / examples section in the appendix (also partly because we don't yet have a proof of correctness, but this is not so critical).
I think it is important to mention this in the main text because NP-DHMC does not work in practice without this adaption.

§ CORRECTNESS

The NP-HMC algorithm is correct in the sense that
the generated Markov chain converges to the target distribution
$\tdist:{A}\mapsto{\frac{1}{Z} \shortint{A}{w}{\tmeasure}}$ [fz]with normalising constant $Z$.
We present an outline of our proof here.
The full proof can be found in
the SM.

Invariant distribution
By iterating <ref>, the NP-HMC algorithm generates a Markov chain $\set{\vec{q}^{(i)}}_{i\in\Nat} $ on the target space $\traces$.
The first step to correctness is to show that the invariant distribution of this chain is indeed the target distribution.

This proof is non-trivial, since the length of the generated sample depends on the values of the random samples drawn in the $\extend$ subroutine (<ref>).
To work around this, we define an auxiliary algorithm,
[lo]which induces the same Markov chain as NP-HMC,
but does not increase the dimension dynamically.
Instead, it finds the smallest $N$ such that all intermediate positions in the $L$ leapfrog steps stay in the domain of $U_N$,
and performs leapfrog steps as in standard HMC.
In this algorithm, all stochastic primitives are executed outside of the Hamiltonian dynamics simulation, and the simulation has a fixed dimension.
Hence we can proceed to identify the invariant distribution.
We then show (in
the SM
that the Markov chain generated by this auxiliary algorithm has the target distribution $\tdist$ as its invariant distribution,
and hence the same holds for NP-HMC (<ref>).

Given Assumptions ass:11, ass:22 and ass:33, the target distribution $\tdist$ is the invariant distribution of the Markov chain generated by iterating <ref>.

the SM
, we extend the proof of <cit.>
to show that the chain converges for a small enough step size $\epsilon$,
as long as the following additional assumptions are met:

* $w$ is continuously differentiable on a non-null set $A$ with measure-zero boundary.

* [fz]$w|_{\support{w}}$ is bounded below by a positive constant.

* For each $n$, the function $\frac{\grad{w_{\le n}}}{w_{\le n}}$ is uniformly bounded from above and below on $\support{w_{\le n}} \cap A$.

* For each $n$, the function $\frac{\grad{w_{\le n}}}{w_{\le n}}$ is Lipschitz [fz]continuous on $\support{w_{\le n}} \cap A$.

If Assumptions (C1)–(C4) are satisfied in addition to Assumptions ass:11, ass:22 and ass:33,
the Markov chain generated by iterating <ref> converges to the target distribution $\tdist$.

@Carol: I tried to phrase the conditions in terms of $w$ instead of $U$. Could you please check this?
Thanks for this.
I fixed condition 2. $w$ can only be positive on the support, so we have to restrict it.

§ EXPERIMENTS

We implemented the NP-HMC algorithm and its variants (NP-RHMC and NP-DHMC) in Python, using PyTorch <cit.> for automatic differentiation.
[fz]We implemented it from scratch rather than in an existing system because NP-DHMC needs additional information about each sample (does the density function depend discontinuously on it?), so it requires a deeper integration in the probabilistic programming system.

In our empirical evaluation, we focus on the NP-DHMC algorithm because it inherits discontinuous HMC's efficient handling of discontinuities:
contrary to NP-RHMC, it does not need to compute the intersections of the particle's trajectory with the regions of discontinuity.
Our implementation also uses the efficiency improvements discussed in
the SM
The code for our implementation and experiments is available at <https://github.com/fzaiser/nonparametric-hmc> and archived as <cit.>.

We compare our implementation with Anglican's <cit.> inference algorithms that are applicable out-of-the-box to nonparametric models:
lightweight Metropolis-Hastings (LMH), particle Gibbs (PGibbs) and random-walk lightweight Metropolis-Hastings (RMH).

[cm]We also compared Interacting Particle MCMC (IPMCMC), but it performed consistently worse than PGibbs in our experiments
[cm]and hence is omitted.
NP-DHMC performs more computation per sample than its competitors because it evaluates the density function in each of the $L$ leapfrog steps, not just once like the other inference algorithms.
To equalise the computation budgets, we generate $L$ times as many samples for each competitor algorithm, and apply thinning (taking every $L$-th sample) to get a comparable sample size.

Total variation distance from the ground truth for the geometric distribution, averaged over 10 runs. Each run: $10^3$ NP-DHMC samples with $10^2$ burn-in, 5 leapfrog steps of size 0.1; and $5\times 10^3$ LMH, PGibbs and RMH samples.

method    ours    LMH    PGibbs    RMH

TVD    0.0136    0.0224    0.0158     0.0196

Geometric distribution
A classic example to illustrate recursion in a universal PPL is sampling from a geometric distribution with parameter $p$ by repeatedly flipping a biased coin with probability $p$ (see e.g. Ch. 5 in <cit.>).
The pseudocode for it is:

def geometric():
if sample(uniform(0, 1)) < p: return 1
else: return 1 + geometric()

We tested our algorithm on this problem with $p = 0.2$.
Our implementation works well on this example and has no trouble jumping between traces of different length.
To quantify this, we computed the total variation distance to the ground truth for each approach and report it in <ref>.
The result is perhaps surprising given that the odds are “stacked against” NP-DHMC in this model:
it is a discrete model, so there is no gradient information, and there are no observations (only sampling from the prior).
These properties should favour the competitors, making the performance of NP-DHMC rather remarkable.

Random walk
To better evaluate our algorithm on probabilistic programs with unbounded loops (such as the example from <ref>), we considered the one-sided random walk on $\pReal$ described in <cit.>:
A pedestrian starts from a random point in $[0,3]$ and walks a uniformly random distance of at most 1 in either direction, until they pass 0.
Given a (noisily) measured total distance of 1.1 travelled, what is the posterior distribution of the starting point?
As this process has infinite expected running time, we need a stopping condition if the pedestrian is too far away from zero, (distance < 10), as shown in the following pseudocode:

start = sample(uniform(0,3))
position = start; distance = 0
while position > 0 and distance < 10:
step = sample(uniform(-1, 1))
position += step; distance += abs(step)
observe(distance, normal(1.1, 0.1))
return start

method     ours    LMH    PGibbs    RMH

ESS    679    526    310     508

Kernel density estimate (top) averaged over 10 runs and estimated effective sample size (bottom) averaged over 10 runs. Each run: $10^3$ NP-DHMC samples with $10^2$ burn-in, 50 leapfrog steps of size 0.1; and $5\times 10^4$ LMH, PGibbs and RMH samples.

This example is interesting and challenging because the true posterior is difficult to determine precisely.
Therefore we took $10^7$ importance samples (effective sample size $\approx 4.4 \times 10^5$) and considered those as the ground truth.
As one can see from <ref>, NP-DHMC comes closest.
Since it is not clear what measure to use for the distance from these “ground truth” samples, we instead computed the effective sample size[We used the standard ESS estimator (based on weighted samples) for the ground-truth importance samples, and an autocorrelation-based MCMC ESS estimator (Sec. 11.5 in <cit.>) for the rest.

for each method (<ref>).
Our method does best in that regard as well.

Kernel density estimate for the random walk example compared to Pyro, averaged over 10 runs. Each run: $10^3$ samples with $10^2$ burn-in, 50 leapfrog steps of size 0.1.

[lo]The popular PPL Pyro accepts nonparametric models as input.

We therefore tried to ascertain the performance of its HMC implementation on this example.
We ran both Pyro's HMC sampler (with the same hyperparameters as ours) and Pyro's No-U-Turn sampler (NUTS) <cit.>, which aims to automatically infer good hyperparameter settings.

The inferred posterior distributions (<ref>) are far away from the ground truth, and clearly wrong.
Pyro’s inference was very slow, sometimes taking almost a minute to produce a single sample.
[lo]For this reason, we didn't run more experiments with it.

Let me know what you think about this. Also, does the plot take up too much space?
I wonder if we could add the Pyro plots to Fig 5. Don't worry if it's too crowded. Another note about colours. I think it would be nice to use a different colour for the Pyro engines and keep the colour of ground truth consistent with Fig 5.
I can do it, but I think it would be too crowded.

Histogram of the no. of mixtures for the GMM example; correct posterior = 9, averaged over 10 runs. Each run: $10^3$ NP-DHMC samples with $10^2$ burn-in, 50 leapfrog steps of size 0.05; and $5\times 10^4$ LMH, PGibbs and RMH samples.

LPPD for the GMM example, averaged over 10 runs. Each run: $10^3$ NP-DHMC samples with $10^2$ burn-in, 50 leapfrog steps of size 0.05; and $5\times 10^4$ LMH, PGibbs and RMH samples. The shaded area is one standard deviation. PGibbs (with final LPPD $-716.85 \pm 0.64$) is omitted to show the top contenders clearly.

LPPD for the DPMM example, averaged over 10 runs. Each run: $10^2$ NP-DHMC samples with $50$ burn-in, 20 leapfrog steps of size 0.05; and $2\times 10^3$ LMH, PGibbs and RMH samples. The shaded area is one standard deviation. PGibbs (with final LPPD $-725.96 \pm 9.83$) is omitted to show the top contenders clearly.

Log pointwise predictive density (LPPD) for the Dirichlet Process Mixture Model (DPMM), mean and standard deviation, averaged over 10 runs. Each run: $10^2$ NP-DHMC samples with $50$ burn-in, 20 leapfrog steps of size 0.05; and $2\times 10^3$ LMH, PGibbs and RMH samples.

model     “true”     ours    LMH    PGibbs    RMH

DPMM    -674.81    $\mathbf{-676.232 \pm 0.62}$    $-679.24 \pm 1.97$    $-725.96 \pm 9.83$    $-679.68 \pm 2.97$

Gaussian mixture model
We also considered the following mixture model adapted from <cit.>, where the number of mixture components $K$ is unbounded.

\begin{align*}
K &\sim \mathrm{Poisson}(10) + 1 \\
\mu_k &\sim \mathrm{Uniform}([0, 100]^3) &&\text{for } k = 1,\dots,K \\
x_n &\sim \frac{1}{K}\sum_{k=1}^K \Gau_3(\mu_k, 10^2I_3) &&\text{for } n = 1,\dots,N
\end{align*}

This model samples parameters $\theta = (K\in \Nat, \mu_k \in [0,100]^3)$ and $N$ data points $X = \{x_1,\dots,x_N\} \subseteq \Real^3$.
Note that this model uses much higher standard deviations (10 instead of 0.1) for the Gaussian mixture components compared to <cit.>.
This is to avoid typical problems of MCMC algorithms with multimodal distributions, which is an issue inherent to MCMC algorithms and tangential to this work.
We used this model to generate $N = 200$ training data points for a fixed $\theta^* = (K^* = 9, \mu_{1\dots K^*}^*)$.
We let our inference algorithms sample from $p(\theta \mid X)$ and compared the posterior on the number of mixture components $K$ with the other inference algorithms.
The histogram (<ref>) shows that NP-DHMC usually finds the correct number of mixture components ($K^* = 9$).

In addition, we computed the log pointwise predictive density (LPPD) for a test set with $N' = 50$ data points $Y = \{y_1,\dots,y_{N'}\}$, generated from the same $\theta^*$ as the training data.
The LPPD is defined as $\sum_{i=1}^{N'} \log \int p(y_i \mid \theta) p(\theta \mid X) \dif\theta$ and can be approximated by $\sum_{i=1}^{N'} \log \frac{1}{M} \sum_{j=1}^M p(y_i \mid \theta_j)$ where $(\theta_j)_{j=1\dots M}$ are samples from $p(\theta \mid X)$ <cit.>.
The results (<ref>) include the “true” LPPD of the test data under the point estimate $\theta^*$.
As we can see, NP-DHMC outperforms the other methods and has the lowest variance over multiple runs.

Dirichlet process mixture model
Finally, we consider a classic example of nonparametric models: the Dirichlet process $\mathrm{DP}(\alpha, H)$ <cit.>, which is a stochastic process parametrised by concentration parameter $\alpha > 0$ and a probability distribution $H$.
For practical purposes, one can think of samples from $\mathrm{DP}(\alpha, H)$ as an infinite sequence $(w_k, h_k)_{n\in \Nat}$ where $w_n$ are weights that sum to 1 and $h_n$ are samples from $H$.
Conceptually, a DP Gaussian mixture model takes the form

\begin{align*}
(w_k, h_k)_{k\in \Nat} &\sim \mathrm{DP}(\alpha, H) \\
x_n &\sim \sum_{k=1}^\infty w_k \cdot \Gau(h_k, \Sigma) &&\text{for } n = 1,\dots,N
\end{align*}

Sampling from a DP is usually implemented by the stick-breaking method <cit.>.
However, one cannot actually compute an infinite sequence.
A practical workaround is to cap the number of mixture components by a fixed $K$ and only consider $(w_k, h_k)_{k=1,\dots,K}$.
[For instance, see <https://pyro.ai/examples/dirichlet_process_mixture.html> (accessed: 2021-06-06), the Pyro tutorial on DP mixture models.]

[lo]This renders the model parametric, enabling the use of standard HMC.
However such a treatment is clearly unsatisfactory: if the data actually requires more mixture components than $K$, the model would be found wanting.

We propose a different approach.
Instead of choosing a fixed $K$, we allow it to depend on the weights $w_k$:
we pick the minimal $K \in \Nat$ such that $\sum_{i=1}^K w_i > 1 - \epsilon$ for some $\epsilon > 0$.
With this restriction, we only discard insignificant mixture components (with a weight $w_k < \epsilon$) and allow as many mixture components as necessary to model the data accurately.
This model is not parametric anymore, but still tractable by our algorithm.

We implemented the above model with the parameters $\alpha = 5$, $\epsilon = 0.01$, and the remaining parameters as chosen in the previous GMM example, i.e. $H = \mathrm{Uniform}([0,100]^3)$ and $\Sigma = 10^2I_3$.
We used the same training and test data as in the previous GMM example and present the LPPD results in <ref>.

As we can see, NP-DHMC outperforms the other methods and has the lowest variance over multiple runs.

§ RELATED WORK AND CONCLUSION

The standard MCMC algorithm for PPLs that is widely implemented (for example, in Anglican, Venture, and Web PPL) is the Lightweight Metropolis-Hastings (LMH) algorithm and its extensions <cit.>,

which performs single-site updates on the current sample and re-executes the program.

Unlike NP-HMC, where Hamiltonian motion is simulated on the resulting extended trace,

LMH suffers from a lack of predictive accuracy in its proposal

(as shown in <ref>).

Due to space constraints, I don't think we need to talk about the MH algorithm from <cit.>.
A lesser-known general purpose inference algorithm for PPLs is the version of Metropolis-Hastings introduced in <cit.>.

In contrast to single-site updates like LMH, it updates all values in the current sample and re-executes the program.

This is the worst performing algorithm in our experiments (see <ref>).

Elementary Inferences This paragraph is standard, and can safely be omitted when space is at a premium.
Inference algorithms like importance sampling (IS) or rejection sampling can be applied to probabilistic programs, and are implemented in most PPLs.
However, as shown in the GMM example in <ref>, IS struggles with high dimension models and
sampling from the prior for the proposal does not give a high acceptance ratio.
Hence, it is not a fitting inference algorithm for complex models described by probabilistic programs.
Todo: Add SMC here

[lo]The Reversible Jump Markov chain Monte Carlo (RJMCMC) algorithm <cit.> is similar to NP-HMC in that it is a trans-dimensional MCMC sampler.

However, NP-HMC is a general purpose inference algorithm that works out-of-the-box when given an input density function, whereas RJMCMC additionally requires the user to specify a transition kernel.
Various RJMCMC transition kernels have been suggested for specific models, e.g. split-merge proposal for infinite Gaussian mixture models.

Some PPLs such as Hakaru, Pyro and Gen give users the flexibility to hand-code the proposal in a MCMC setting.
For instance, <cit.> implement the split-merge proposal <cit.> of RJMCMC in Gen.
Though this line of research is orthogonal to ours,
PPLs such as Gen could play a useful role in the implementation of NP-HMC and similar extensions of inference algorithms to nonparametric models.
This is however not ideal as designing efficient proposals for specific complex probabilistic models requires knowledge in the workings of inference algorithms and is usually not straightforward.
Moreover, this setting forfeits one of the main appeals of PPLs to automate Bayesian inference.
We should be a little more positive in our assessment. This line of work is well-motivated (but orthogonal to ours); it is in our interest to acknowledge it.

The HMC algorithm and its variants, notably the No-U-Turn Sampler,

are the workhorse inference methods in the influential PPL Stan <cit.>.

The challenges posed by stochastic branching in PPLs are the focus of reflective/refractive HMC <cit.>; discontinuous HMC <cit.>; mixed HMC <cit.>; and the first-order PPL in <cit.> which is equipped with an implementation of discontinuous HMC.
[lo]By contrast, our work is an attempt to tackle the language constructs of branching and recursion.

[lo]Unlike Monte Carlo methods, variational inference (VI) <cit.> solves the Bayesian inference problem by treating it as an optimisation problem.
When adapted to models expressed as probabilistic programs, the score function VI <cit.> can in principle be applied to a large class of branching and recursive programs because only the variational density functions need to be differentiable.
Existing implementations of VI algorithms in probabilistic programming systems are however far from automatic:
in the main, the guide programs (that express variational distributions) still need to be hand-coded.

A different approach to posterior inference is Stochastic Variational Inference (SVI) <cit.>:
the goal is to approximate the posterior distribution with a parametrised family of simpler distributions.
Thus, one wants to find the parameters that are closest to the posterior according to some distance measure (e.g. Kullback–Leibler divergence).
This optimization problem can be approximately solved using stochastic gradient descent.

[cm]We think SVI and NP-HMC would complement each other in PPLs with SVI more fitting for models with large dataset and NP-HMC superior when the dataset is small.
Should we just say that MCMC and SVI are used in different settings? MCMC is better when the dataset is small and accuracy is required whereas SVI is better for a large dataset.

Recently, <cit.> introduced the Divide, Conquer, and Combine (DCC) algorithm, which is applicable to programs definable using branching and recursion.
[fz]As a hybrid algorithm, DCC solves the problem of designing a proposal that can efficiently transition between configurations by performing local inferences on submodels, and returning an appropriately weighted combination of the respective samples.
Thanks to a judicious resource allocation scheme, it exhibits strong performance on multimodal distributions.

Recently, <cit.> introduced the Divide, Conquer, and Combine (DCC) algorithm, which performs local inferences on submodels of an input probabilistic program, and returns a combination of the respective samples.

[cm]DCC can be seen as a hybrid method which our NP-HMC algorithm might complement in a universal PPL.

Do they say this in the DCC paper? We should only claim this if we know for sure.
Thanks for checking my worst instincts. I always feel like I need to address the shortcomings of other algorithms.
Zhou: This is not fully true though, as LMH forgets any proposed, not accepted samples whereas DCC ‘remembers’ them. Also DCC uses RMH actually.

I agree that it is possible to make proposals cleverer. However, at the end of the day, if we are not provided more information, we will need to rely on sth. similar to RMH to explore (prior + local movements).

I cannot comment for now as I do not fully follow how NP-HMC extends. But I think it might be good to compare your jumping scheme with LMH, in order to improve clarity.

The paper that inspired this work most is the probabilistic path HMC <cit.> which performs HMC on discrete tree spaces described by orthant complexes.
Though addressing a very different problem, its treatment of movements between tree topologies helped us to understand movement between branches in the program tree of a general probabilistic program.
This is ultimately not so relevant. It can go if space is an issue.

§.§ Conclusion

[lo]We have presented the NP-HMC algorithm, the first extension of the HMC algorithm to nonparametric models.

We have proved that NP-HMC is correct. We have also empirically demonstrated that it enjoys significant performance improvements over state-of-the-art MCMC algorithms for universal PPLs on four nonparametric models, thereby illustrating that the key advantage of HMC—the proposal of moves to distant states with a high acceptance probability—has been preserved by NP-HMC.

[fz]We have presented the nonparametric Hamiltonian Monte Carlo algorithm, which generalises the popular HMC algorithm to nonparametric models typically described by universal PPLs.
We highlight its theoretical and practical merits in preserving the effectiveness of HMC and outperforming state-of-the-art MCMC algorithms for universal PPLs.

We have presented an extension of the popular HMC algorithm to nonparametric models.
NP-HMC enjoys nice theoretical properties: the adaption from standard HMC is conceptually simple (one only has to insert $\extend$ calls in the integrator), and we have proved its correctness.
In practice, our algorithm inherits the effectiveness of HMC and outperforms state-of-the-art inference algorithms for universal PPLs.

We thank the reviewers for their insightful feedback and pointing out [fz]important related work.
We are grateful to Hugo Paquet, Dominik Wagner and Yuan Zhou, who gave detailed comments on an early draft, and to Tom Rainforth and Arnaud Doucet for their helpful comments and advice.
[fz]We gratefully acknowledge support from the EPSRCNot sure how to word this, but EPSRC is funding my DPhil.[cm] and the Croucher Foundation.
Thanks to funding agencies and corporate sponsors that provided financial support?

Nonparametric Hamiltonian Monte Carlo (Appendix)

§ STATISTICAL PCF

§ HAMILTONIAN MONTE CARLO ALGORITHM AND ITS VARIANTS

§ PROOF OF CORRECTNESS

§ EXPERIMENTS

§ DO NOT HAVE AN APPENDIX HERE

Do not put content after the references.

Put anything that you might normally include after the references in a separate
supplementary file.

We recommend that you build supplementary material in a separate document.
If you must create one PDF and cut it up, please be careful to use a tool that
doesn't alter the margins, and that doesn't aggressively rewrite the PDF file.
pdftk usually works fine.

Please do not use Apple's preview to cut off supplementary material. In
previous years it has altered margins, and created headaches at the camera-ready