2017; Sereno et al., 2017; McClintock et al., 2019). Sources of systematics
can be taken as uncorrelated or correlated in different analyses (Melchior et
al., 2017). Notwithstanding the differences, different analyses agree on the
main sources of systematics and on their size. Mass biases of the order of
5–10 % are consistent with results from our cross-comparison.

We perform an end-to-end analysis, and we compare the main results, which are
affected by all potential systematics, which enter at different steps.
However, cross-checks can be also used to assess specific systematic effects.
In App. C, we show how residual effects in shear calibration for Stage-III
surveys and precursors are negligible with respect to the statistical
uncertainties by comparing the WL profiles of clusters extracted from matched
source catalogues.

##  12 Conclusions

Both statistical uncertainty and unaccounted for systematics can limit the
accuracy of the mass calibration. The main contributors to the systematic
error budget are the calibration uncertainties of the shear measurements, the
photo-$z$ performance, and the selection of the source galaxies. Furthermore,
in the regime of cluster lensing, shape calibration has to be revised for the
impact of stronger shears and increased blending (Hernández-Martín et al.,
2020). Köhlinger et al. (2015) explored statistical uncertainties and
systematic errors in WL mass estimates of galaxy clusters in a Euclid-like
survey and showed that requirements for WL by clusters are similar but a bit
weaker than what is needed for cosmic shear. In principle, proper calibration
for cosmic shear should also meet requirements for WL clusters.

To test the residual level of systematic effects in Stage-III surveys, we
performed a cross-survey analysis. To assess the robustness of the cluster
mass determinations, we performed a unified analysis exploiting the COMB-CL
processing function to be used by the Euclid Collaboration for WL cluster mass
measurements. On one hand, COMB-CL can deal with different data sets, or shear
and photo-$z$ estimation algorithms, and, on the other hand, it can measure
the WL signal and determine the cluster mass with uniform and consistent model
assumptions, and the same inference pipeline. The analysis also served as a
testing of the COMB-CL processing function.

Our choices and assumptions were meant to provide a consistent framework for
comparison of very different data sets. The goal was to look for
discrepancies, not to rank alternative methodologies for WL cluster analyses
or surveys, or fully exploit their statistical power. For example, we adopted
the same secure colour-colour cut for all surveys. This made the background
selection and the systematic effects related to member and foreground
contamination quite consistent. Different cuts which were optimised per survey
could improve the purity and completeness of background selection at a price
of a less straightforward cross-survey comparison.

The comparison of the cluster lensing signal and of the inferred WL mass
showed a remarkable agreement between independent imaging surveys. Differences
in excess surface density profiles were consistent with the statistical
uncertainties and well within survey defined systematic error budgets for
Stage-III surveys and precursors. This is consistent with previous results on
cosmic shear (Chang et al., 2019; Longley et al., 2023), or galaxy-galaxy
lensing (Leauthaud et al., 2022). Masses of individual clusters were found to
agree with literature results. Our analysis suggests that unknown or residual
systematic effects are subdominant with respect to statistical noise or known
systematic errors for Stage-III surveys and precursors.

However, in order to improve WL mass calibration, difficult challenges remain.
Shear calibration can be hampered by incomplete image simulations (Mandelbaum
et al., 2018). Sufficient archival Hubble Space Telescope data should be
available to create realistic populations of mock galaxies for Euclid-like
simulations and achieve the required sub-percent level accuracy (Hoekstra et
al., 2017). Performance could also improve with data-driven approaches such as
metacalibration (Huff & Mandelbaum, 2017; Sheldon & Huff, 2017; Hoekstra,
2021).

Spectroscopic training samples used to calibrate photo-$z$s can be incomplete
at high redshifts or for some galaxy types, resulting in systematic
uncertainty in the photo-$z$ calibration. With these issues in mind, the
Euclid Collaboration has been making a significant effort to enlarge the
calibration sample (see, e.g., Euclid Collaboration: Saglia et al., 2022) and
properly train the photo-$z$ algorithms (see, e.g., Euclid Collaboration et
al., 2023).

In any analysis, one seeks to obtain the right balance between high
statistical precision and a low level of systematic errors. For example, we
favoured very conservative background galaxy selections, where we retained
only the sources with high fidelity redshift estimates. On one hand, this
results in a larger shape noise than in, e.g., analyses using the full source
sample. On the other hand, systematic uncertainties are expected to be
smaller.

Precise mass calibration is required to ensure the success of Stage-IV cluster
cosmology measurements (Sartoris et al., 2016). Mass calibration at the level
of a few percent is already within reach of Stage-III analyses. In a Euclid-
like survey, shear is going to be calibrated within an accuracy of $|\Delta
m|\sim 2\times 10^{-3}$ (Cropper et al., 2013). Given the exquisite space-
based infrared photometry complemented by ground-based optical coverage and
spectroscopic calibration samples (Euclid Collaboration: Saglia et al., 2022),
photo-$z$ estimates will be very accurate, with an expected bias of the order
of $|\Delta z_{\text{p}}|\sim 0.002(1+z)$ (Euclid Collaboration et al., 2023).
Such accurate photo-$z$s along with extended optical/infrared photometry will
also strongly aid background selection. We expect to select nearly complete
(up to $\sim 90\,\%$) and pure samples of background sources, with a
contamination at the sub-percent level (Euclid Collaboration: Lesci et al.,
2023). Systematics due to analysis assumptions can be significantly reduced
with lens modelling accounting for, e.g., miscentring, concentration and
truncation, and the effect of correlated or uncorrelated matter (Sereno et
al., 2018; Bellagamba et al., 2019; McClintock et al., 2019). Mass calibration
at the percent level or better could be within the reach of Stage-IV surveys.

###### Acknowledgements.

This work was made possible by public sharing of tools and data sets by the
astronomical community, which we thanks in the following with the suggested
formulations. MS acknowledges financial contributions from contract ASI-INAF
n.2017-14-H.0, contract INAF mainstream project 1.05.01.86.10, INAF Theory
Grant 2023: Gravitational lensing detection of matter distribution at galaxy
cluster boundaries and beyond (1.05.23.06.17), and contract Prin-MUR 2022
supported by Next Generation EU (n.20227RNLY3, The concordance cosmological
model: stress-tests with galaxy clusters). CG thanks the support from INAF
theory Grant 2022: Illuminating Dark Matter using Weak Lensing by Cluster
Satellites, PI: Carlo Giocoli. LM acknowledges support from the grants PRIN-
MIUR 2017 WSCC32 and ASI n.2018-23-HH.0. The Euclid Consortium acknowledges
the European Space Agency and a number of agencies and institutes that have
supported the development of Euclid, in particular the Academy of Finland, the
Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid
Consortium, the French Centre National d’Etudes Spatiales, the Deutsches
Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the
Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e
Innovación, the National Aeronautics and Space Administration, the National
Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor
Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State
Secretariat for Education, Research and Innovation (SERI) at the Swiss Space
Office (SSO), and the United Kingdom Space Agency. A complete and detailed
list is available on the Euclid web site (http://www.euclid-ec.org). This
research has made use of NASA’s Astrophysics Data System (ADS) and of the
NASA/IPAC Extragalactic Database (NED), which is operated by the Jet
Propulsion Laboratory, California Institute of Technology, under contract with
the National Aeronautics and Space Administration. The Hyper Suprime-Cam (HSC)
collaboration includes the astronomical communities of Japan and Taiwan, and
Princeton University. The HSC instrumentation and software were developed by
the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for
the Physics and Mathematics of the Universe (Kavli IPMU), the University of
Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia
Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and
Princeton University. Funding was contributed by the FIRST program from
Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science
and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS),
Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ,
Kavli IPMU, KEK, ASIAA, and Princeton University. Based in part on data
collected at the Subaru Telescope and retrieved from the HSC data archive
system, which is operated by Subaru Telescope and Astronomy Data Center at
National Astronomical Observatory of Japan. This work is based in part on
observations obtained with MegaPrime/MegaCam, a joint project of CFHT and
CEA/IRFU, at the Canada-France-Hawaii Telescope (CFHT) which is operated by
the National Research Council (NRC) of Canada, the Institut National des
Sciences de l’Univers of the Centre National de la Recherche Scientifique
(CNRS) of France, and the University of Hawaii. This research used the
facilities of the Canadian Astronomy Data Centre operated by the National
Research Council of Canada with the support of the Canadian Space Agency.
CFHTLenS and RCSLenS data processing was made possible thanks to significant
computing support from the NSERC Research Tools and Instruments grant program.
This work is based in part on observations made with ESO Telescopes at the La
Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017,
177.A-3018 and 179.A-2004, and on data products produced by the KiDS
consortium. The KiDS production team acknowledges support from: Deutsche
Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of
Padova, and the University Federico II (Naples). This project used public
archival data from the Dark Energy Survey (DES). Funding for the DES Projects
has been provided by the U.S. Department of Energy, the U.S. National Science
Foundation, the Ministry of Science and Education of Spain, the Science and
Technology FacilitiesCouncil of the United Kingdom, the Higher Education
Funding Council for England, the National Center for Supercomputing
Applications at the University of Illinois at Urbana-Champaign, the Kavli
Institute of Cosmological Physics at the University of Chicago, the Center for
Cosmology and Astro-Particle Physics at the Ohio State University, the
Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M
University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho
de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de
Desenvolvimento Científico e Tecnológico and the Ministério da Ciência,
Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the
Collaborating Institutions in the Dark Energy Survey. The Collaborating
Institutions are Argonne National Laboratory, the University of California at
Santa Cruz, the University of Cambridge, Centro de Investigaciones
Energéticas, Medioambientales y Tecnológicas-Madrid, the University of
Chicago, University College London, the DES-Brazil Consortium, the University
of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi
National Accelerator Laboratory, the University of Illinois at Urbana-
Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de
Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-
Maximilians Universität München and the associated Excellence Cluster
Universe, the University of Michigan, the National Optical Astronomy
Observatory, the University of Nottingham, The Ohio State University, the
OzDES Membership Consortium, the University of Pennsylvania, the University of
Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the
University of Sussex, and Texas A&M University. Based in part on observations
at Cerro Tololo Inter-American Observatory, National Optical Astronomy
Observatory, which is operated by the Association of Universities for Research
in Astronomy (AURA) under a cooperative agreement with the National Science
Foundation.

## References

  * Abbott et al. (2020) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, Phys. Rev. D, 102, 023509
  * Adami et al. (2018) Adami, C., Giles, P., Koulouridis, E., et al. 2018, A&A, 620, A5
  * Aihara et al. (2018) Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S8
  * Albrecht et al. (2006) Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, arXiv:astro-ph/0609591
  * Andreon & Bergé (2012) Andreon, S. & Bergé, J. 2012, A&A, 547, A117
  * Applegate et al. (2014) Applegate, D. E., von der Linden, A., Kelly, P. L., et al. 2014, MNRAS, 439, 48
  * Asgari et al. (2021) Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104
  * Baltz et al. (2009) Baltz, E. A., Marshall, P., & Oguri, M. 2009, J. Cosmology Astropart. Phys., 1, 15
  * Bartelmann & Schneider (2001) Bartelmann, M. & Schneider, P. 2001, Phys. Rep, 340, 291
  * Beers et al. (1990) Beers, T. C., Flynn, K., & Gebhardt, K. 1990, AJ, 100, 32
  * Bellagamba et al. (2018) Bellagamba, F., Roncarelli, M., Maturi, M., & Moscardini, L. 2018, MNRAS, 473, 5221
  * Bellagamba et al. (2019) Bellagamba, F., Sereno, M., Roncarelli, M., et al. 2019, MNRAS, 484, 1598
  * Benítez (2000) Benítez, N. 2000, ApJ, 536, 571
  * Benjamin et al. (2013) Benjamin, J., Van Waerbeke, L., Heymans, C., et al. 2013, MNRAS, 431, 1547
  * Bernstein & Jarvis (2002) Bernstein, G. M. & Jarvis, M. 2002, AJ, 123, 583
  * Bocquet et al. (2019) Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55
  * Chang et al. (2019) Chang, C., Wang, M., Dodelson, S., et al. 2019, MNRAS, 482, 3696
  * CHEX-MATE Collaboration et al. (2021) CHEX-MATE Collaboration, Arnaud, M., Ettori, S., et al. 2021, A&A, 650, A104
  * Chiu et al. (2020) Chiu, I. N., Umetsu, K., Murata, R., Medezinski, E., & Oguri, M. 2020, MNRAS, 495, 428
  * Costanzi et al. (2021) Costanzi, M., Saro, A., Bocquet, S., et al. 2021, Phys. Rev. D, 103, 043522
  * Cropper et al. (2013) Cropper, M., Hoekstra, H., Kitching, T., et al. 2013, MNRAS, 431, 3103
  * Dark Energy Survey et al. (2023) Dark Energy Survey, Kilo-Degree Survey Collaboration, Abbott, T. M. C., et al. 2023, arXiv:2305.17173
  * de Jong et al. (2013) de Jong, J. T. A., Kuijken, K., Applegate, D., et al. 2013, The Messenger, 154, 44
  * Despali et al. (2016) Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486
  * Diemer & Joyce (2019) Diemer, B. & Joyce, M. 2019, ApJ, 871, 168
  * Diemer & Kravtsov (2014) Diemer, B. & Kravtsov, A. V. 2014, ApJ, 789, 1
  * Dietrich et al. (2019) Dietrich, J. P., Bocquet, S., Schrabback, T., et al. 2019, MNRAS, 483, 2871
  * Eisenstein & Hu (1998) Eisenstein, D. J. & Hu, W. 1998, ApJ, 496, 605
  * Erben et al. (2013) Erben, T., Hildebrandt, H., Miller, L., et al. 2013, MNRAS, 433, 2545
  * Ettori et al. (2009) Ettori, S., Morandi, A., Tozzi, P., et al. 2009, A&A, 501, 61
  * Euclid Collaboration et al. (2023) Euclid Collaboration, Bisigello, L., Conselice, C. J., et al. 2023, MNRAS, 520, 3529
  * Euclid Collaboration: Adam et al. (2019) Euclid Collaboration: Adam, R., Vannier, M., Maurogordato, S., et al. 2019, A&A, 627, A23
  * Euclid Collaboration: Ajani et al. (2023) Euclid Collaboration: Ajani, V., Baldi, M., Barthelemy, A., et al. 2023, A&A, 675, A120
  * Euclid Collaboration: Lesci et al. (2023) Euclid Collaboration: Lesci, G. F., Sereno, M., Radovich, M., et al. 2023, arXiv:2311.16239
  * Euclid Collaboration: Martinet et al. (2019) Euclid Collaboration: Martinet, N., Schrabback, T., Hoekstra, H., et al. 2019, A&A, 627, A59
  * Euclid Collaboration: Saglia et al. (2022) Euclid Collaboration: Saglia, R., De Nicola, S., Fabricius, M., et al. 2022, A&A, 664, A196
  * Euclid Collaboration: Scaramella et al. (2022) Euclid Collaboration: Scaramella, R., Amiaux, J., Mellier, Y., et al. 2022, A&A, 662, A112
  * Fenech Conti et al. (2017) Fenech Conti, I., Herbonnet, R., Hoekstra, H., et al. 2017, MNRAS, 467, 1627
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fu et al. (2008) Fu, L., Semboloni, E., Hoekstra, H., et al. 2008, A&A, 479, 9
  * Furusawa et al. (2018) Furusawa, H., Koike, M., Takata, T., et al. 2018, PASJ, 70, S3
  * Gatti et al. (2021) Gatti, M., Sheldon, E., Amon, A., et al. 2021, MNRAS, 504, 4312
  * Giblin et al. (2021) Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105
  * Gilbank et al. (2011) Gilbank, D. G., Gladders, M. D., Yee, H. K. C., & Hsieh, B. C. 2011, AJ, 141, 94
  * Grandis et al. (2021) Grandis, S., Bocquet, S., Mohr, J. J., Klein, M., & Dolag, K. 2021, MNRAS, 507, 5671
  * Hamana et al. (2020) Hamana, T., Shirasaki, M., Miyazaki, S., et al. 2020, PASJ, 72, 16
  * Hernández-Martín et al. (2020) Hernández-Martín, B., Schrabback, T., Hoekstra, H., et al. 2020, A&A, 640, A117
  * Heymans et al. (2021) Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140
  * Heymans et al. (2012) Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146
  * Hikage et al. (2019) Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43
  * Hildebrandt et al. (2016) Hildebrandt, H., Choi, A., Heymans, C., et al. 2016, MNRAS, 463, 635
  * Hildebrandt et al. (2012) Hildebrandt, H., Erben, T., Kuijken, K., et al. 2012, MNRAS, 421, 2355
  * Hildebrandt et al. (2020) Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69
  * Hildebrandt et al. (2021) Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124
  * Hirata & Seljak (2003) Hirata, C. & Seljak, U. 2003, MNRAS, 343, 459
  * Hoekstra (2003) Hoekstra, H. 2003, MNRAS, 339, 1155
  * Hoekstra (2021) Hoekstra, H. 2021, A&A, 656, A135
  * Hoekstra et al. (2015) Hoekstra, H., Herbonnet, R., Muzzin, A., et al. 2015, MNRAS, 449, 685
  * Hoekstra et al. (2012) Hoekstra, H., Mahdavi, A., Babul, A., & Bildfell, C. 2012, MNRAS, 427, 1298
  * Hoekstra et al. (2017) Hoekstra, H., Viola, M., & Herbonnet, R. 2017, MNRAS, 468, 3295
  * Huang et al. (2018) Huang, S., Leauthaud, A., Murata, R., et al. 2018, PASJ, 70, S6
  * Huff & Mandelbaum (2017) Huff, E. & Mandelbaum, R. 2017, arXiv:1702.02600
  * Jarvis et al. (2016) Jarvis, M., Sheldon, E., Zuntz, J., et al. 2016, MNRAS, 460, 2245
  * Joachimi et al. (2021) Joachimi, B., Lin, C. A., Asgari, M., et al. 2021, A&A, 646, A129
  * Jullo et al. (2010) Jullo, E., Natarajan, P., Kneib, J.-P., et al. 2010, Science, 329, 924
  * Kaiser et al. (1995) Kaiser, N., Squires, G., & Broadhurst, T. 1995, ApJ, 449, 460
  * Kannawadi et al. (2019) Kannawadi, A., Hoekstra, H., Miller, L., et al. 2019, A&A, 624, A92
  * Kawanomoto et al. (2018) Kawanomoto, S., Uraguchi, F., Komiyama, Y., et al. 2018, PASJ, 70, 66
  * Köhlinger et al. (2015) Köhlinger, F., Hoekstra, H., & Eriksen, M. 2015, MNRAS, 453, 3107
  * Komiyama et al. (2018) Komiyama, Y., Obuchi, Y., Nakaya, H., et al. 2018, PASJ, 70, S2
  * Kuijken et al. (2019) Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193
  * Leauthaud et al. (2022) Leauthaud, A., Amon, A., Singh, S., et al. 2022, MNRAS, 510, 6150
  * Lesci et al. (2022) Lesci, G. F., Marulli, F., Moscardini, L., et al. 2022, A&A, 659, A88
  * Longley et al. (2023) Longley, E. P., Chang, C., Walter, C. W., et al. 2023, MNRAS, 520, 5016
  * Lubini et al. (2014) Lubini, M., Sereno, M., Coles, J., Jetzer, P., & Saha, P. 2014, MNRAS, 437, 2461
  * Mandelbaum et al. (2018) Mandelbaum, R., Miyatake, H., Hamana, T., et al. 2018, PASJ, 70, S25
  * Mandelbaum et al. (2008) Mandelbaum, R., Seljak, U., & Hirata, C. M. 2008, J. Cosmology Astropart. Phys., 8, 006
  * Mantz et al. (2010) Mantz, A., Allen, S. W., Rapetti, D., & Ebeling, H. 2010, MNRAS, 406, 1759
  * Martinet et al. (2021) Martinet, N., Harnois-Déraps, J., Jullo, E., & Schneider, P. 2021, A&A, 646, A62
  * Maturi et al. (2019) Maturi, M., Bellagamba, F., Radovich, M., et al. 2019, MNRAS, 485, 498
  * McClintock et al. (2019) McClintock, T., Varga, T. N., Gruen, D., et al. 2019, MNRAS, 482, 1352
  * Medezinski et al. (2018) Medezinski, E., Oguri, M., Nishizawa, A. J., et al. 2018, PASJ, 70, 30
  * Melchior et al. (2017) Melchior, P., Gruen, D., McClintock, T., et al. 2017, MNRAS, 469, 4899
  * Melchior et al. (2015) Melchior, P., Suchyta, E., Huff, E., et al. 2015, MNRAS, 449, 2219
  * Miller et al. (2013) Miller, L., Heymans, C., Kitching, T. D., et al. 2013, MNRAS, 429, 2858
  * Miyatake et al. (2019) Miyatake, H., Battaglia, N., Hilton, M., et al. 2019, ApJ, 875, 63
  * Miyatake et al. (2013) Miyatake, H., Nishizawa, A. J., Takada, M., et al. 2013, MNRAS, 429, 3627
  * Miyazaki et al. (2018) Miyazaki, S., Komiyama, Y., Kawanomoto, S., et al. 2018, PASJ, 70, S1
  * Murray et al. (2022) Murray, C., Bartlett, J. G., Artis, E., & Melin, J.-B. 2022, MNRAS, 512, 4785
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
  * Oguri (2014) Oguri, M. 2014, MNRAS, 444, 147
  * Oguri & Hamana (2011) Oguri, M. & Hamana, T. 2011, MNRAS, 414, 1851
  * Oguri et al. (2018) Oguri, M., Lin, Y.-T., Lin, S.-C., et al. 2018, PASJ, 70, S20
  * Okabe & Smith (2016) Okabe, N. & Smith, G. P. 2016, MNRAS, 461, 3794
  * Pierre et al. (2016) Pierre, M., Pacaud, F., Adami, C., et al. 2016, A&A, 592, A1
  * Planck Collaboration: Ade et al. (2014) Planck Collaboration: Ade, P. A. R., Aghanim, N., Armitage-Caplan, C., et al. 2014, A&A, 571, A20
  * Planck Collaboration: Ade et al. (2016) Planck Collaboration: Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2016, A&A, 594, A27
  * Planck Collaboration: Aghanim et al. (2020) Planck Collaboration: Aghanim, N., Akrami, Y., Ashdown, M., et al. 2020, A&A, 641, A6
  * Rowe et al. (2015) Rowe, B. T. P., Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Computing, 10, 121
  * Rykoff et al. (2014) Rykoff, E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104
  * Rykoff et al. (2016) Rykoff, E. S., Rozo, E., Hollowood, D., et al. 2016, ApJS, 224, 1
  * Sartoris et al. (2016) Sartoris, B., Biviano, A., Fedeli, C., et al. 2016, MNRAS, 459, 1764
  * Schneider et al. (1998) Schneider, P., van Waerbeke, L., Jain, B., & Kruse, G. 1998, MNRAS, 296, 873
  * Schrabback et al. (2018a) Schrabback, T., Applegate, D., Dietrich, J. P., et al. 2018a, MNRAS, 474, 2635
  * Schrabback et al. (2018b) Schrabback, T., Schirmer, M., van der Burg, R. F. J., et al. 2018b, A&A, 610, A85
  * Schramm & Kayser (1995) Schramm, T. & Kayser, R. 1995, A&A, 299, 1
  * Seitz & Schneider (1997) Seitz, C. & Schneider, P. 1997, A&A, 318, 687
  * Sereno (2002) Sereno, M. 2002, A&A, 393, 757
  * Sereno (2015) Sereno, M. 2015, MNRAS, 450, 3665
  * Sereno et al. (2017) Sereno, M., Covone, G., Izzo, L., et al. 2017, MNRAS, 472, 1946
  * Sereno et al. (2016) Sereno, M., Fedeli, C., & Moscardini, L. 2016, J. Cosmology Astropart. Phys., 1, 042
  * Sereno et al. (2018) Sereno, M., Giocoli, C., Izzo, L., et al. 2018, Nature Astronomy, 2, 744
  * Sereno et al. (2015) Sereno, M., Veropalumbo, A., Marulli, F., et al. 2015, MNRAS, 449, 4147
  * Sheldon & Huff (2017) Sheldon, E. S. & Huff, E. M. 2017, ApJ, 841, 24
  * Sheldon et al. (2004) Sheldon, E. S., Johnston, D. E., Frieman, J. A., et al. 2004, AJ, 127, 2544
  * Simet et al. (2017) Simet, M., McClintock, T., Mandelbaum, R., et al. 2017, MNRAS, 466, 3103
  * Takahashi et al. (2012) Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152
  * Tanaka et al. (2018) Tanaka, M., Coupon, J., Hsieh, B.-C., et al. 2018, PASJ, 70, S9
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709
  * Tinker et al. (2010) Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878
  * Tröster et al. (2021) Tröster, T., Asgari, M., Blake, C., et al. 2021, A&A, 649, A88
  * Umetsu (2020) Umetsu, K. 2020, A&A Rev., 28, 7
  * Umetsu et al. (2014) Umetsu, K., Medezinski, E., Nonino, M., et al. 2014, ApJ, 795, 163
  * Umetsu et al. (2020) Umetsu, K., Sereno, M., Lieu, M., et al. 2020, ApJ, 890, 148
  * Umetsu et al. (2016) Umetsu, K., Zitrin, A., Gruen, D., et al. 2016, ApJ, 821, 116
  * Viola et al. (2015) Viola, M., Cacciato, M., Brouwer, M., et al. 2015, MNRAS, 452, 3529
  * Voit (2005) Voit, G. M. 2005, Reviews of Modern Physics, 77, 207
  * von der Linden et al. (2014) von der Linden, A., Allen, M. T., Applegate, D. E., et al. 2014, MNRAS, 439, 2

##  Appendix A COMB-CL Pipeline

Figure 10: Structure of the Euclid combined cluster and weak-lensing pipeline,
COMB-CL. White parallelograms show input or output products and blue boxes
show COMB-CL modules that perform different processing steps. Table 15: Survey
specifics to instantiate the COMB-CL Shear class

option | HSC-SSP S16a | CFHTLenS | KiDS DR4 | DES SV1 | RCSLenS
---|---|---|---|---|---
shape_estimator | distortion | ellipticity | ellipticity | ellipticity | ellipticity
sky_coordinates | eastward | westward | westward | westward | westward
selection_photo_z | True | True | True | True | True
z_p_fails | [-99,99] | [-99,99] | [-99,99] | [-99,99] | [-99,99]
Delta_z_lens_0 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1
Delta_z_lens_z | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
z_p_range_min | 0.2 | 0.2 | 0.2 | 0.2 | 0.4
z_p_range_max | 1.5 | 1.2 | 1.2 | 1.2 | 1.2
z_p_conf_min | 0.8 | 0.8 | 0.8 | 0.8 | 0.8
colour_cut_method | ’Medezinski18_griz’ | ’Medezinski18_griz’ | ’Medezinski18_griz’ | ’Medezinski18_griz’ | ’Medezinski18_griz’
mag_fails | [-99,99] | [-99,99] | [-99,99] | [-99,99] | [-99,99]
alpha | 1.0 | 1.0 | 1.0 | 1.0 | 1.0
delta_alpha | 1.0 | 1.0 | 1.0 | 1.0 | 1.0

The analysis presented in this paper is performed with the combined cluster
and weak-lensing pipeline, COMB-CL, that forms part of the global Euclid data
processing pipeline. In the following text, we provide a brief overview of the
structure of COMB-CL as well as the specific configuration choices made for
this work. A more comprehensive description of the code and methodology
implemented in COMB-CL will be provided with the code public release.

The structure of COMB-CL is presented in Fig. 10, where input and output
products are shown in white parallelograms, while the various modules that
perform different processing steps are shown in blue boxes. The input cluster
catalogue is expected to contain the positions, redshifts and, if available,
richness estimates for the detected objects. The input galaxy catalogue is
expected to contain the positions of source galaxies with measured shear
values, photometric redshift point estimates and, if available, magnitudes or
fluxes.

The Cosmology module simply defines a shared class instance for the
cosmological parameters provided ($H_{0}$, $\Omega_{\text{b}}$,
$\Omega_{\text{m}}$, $n_{\text{s}}$ and $\sigma_{8}$). This instance provides
consistent cosmological parameters and functions, such as the Hubble
parameter, $H(z)$, and the comoving distance, $d_{\text{C}}(z)$, to the other
COMB-CL modules.

The Halo Bias, and Halo Concentration modules implement various models for the
halo bias and concentration as a function of mass, respectively. The models
used for the work presented in this paper are the halo bias model of Tinker et
al. (2010) and the halo concentration model of Diemer & Joyce (2019).

The Halo Module module takes a given combination of halo and concentration
models and defines various models of the 1-halo and 2-halo components of the
surface mass density as a function of radius. These components can then be
combined to provide a final model for the reduced differential density
profile. The specific model used for the work presented in this paper is
described in Sect. 6.1.

The Shear Profile module reads in the galaxy and cluster catalogues provided.
Then, for each cluster, background source galaxies are selected using a series
of cuts in photometric redshift and/or colour as described in Sect. 4. The
tangential shear profile for each cluster is then derived as described in
Sect. 5 and provided as an output product of COMB-CL. In Table 15, we list the
specific argument values passed to the Shear Profile module for this work.

Finally, the MCMC Fitting module performs a Markov chain Monte Carlo (MCMC)
analysis to fit the defined halo model to the measured tangential shear
profiles to derive weak-lensing mass estimates, as well as other parameters
describing the 1-halo term (e.g. concentration) or the 2-halo term (e.g., halo
bias), for each of the clusters. The fitting procedure is described in Sect.
6.2. The specific package employed by COMB-CL to perform the MCMC analysis is
emcee (Foreman-Mackey et al. 2013).

##  Appendix B Planck clusters in multiple surveys

---
Figure 11: Reduced excess surface surface density profiles of Planck clusters
detected in multiple public surveys. Table 16: Planck clusters in multiple
surveys. Comparison of the shear profiles as derived from survey $a$ (col. 2)
or $b$ (col. 3); $\chi^{2}_{a,b}$ is the $\chi^{2}$ between the two profiles;
$p_{a,b}=p\left(\chi^{2}>\chi_{a,b}^{2}\right)$ is the probability to exceed
this value.

Cluster | survey a | survey b | $\chi^{2}_{a,b}$ | $p_{a,b}$
---|---|---|---|---
PSZ2 G167.98$-$59.95 | HSC-SSP S16a | CFHTLenS | 2.3 | 0.97
PSZ2 G174.40$-$57.33 | HSC-SSP S16a | CFHTLenS | 3.0 | 0.93
PSZ2 G230.73$+$27.70 | CFHTLenS | KiDS DR4 | 14.7 | 0.07
PSZ2 G231.79$+$31.48 | HSC-SSP S16a | KiDS DR4 | 8.6 | 0.38

For four PSZ2 clusters, we retrieve data from two surveys. The measured
reduced excess surface density profiles are shown in Fig. 11, whereas the
statistical comparison is summarised in Table 16. We find no evidence for
disagreement.

##  Appendix C Shear calibration

Table 17:  Systematic error budget on the shear calibration for WL cluster
analyses in Stage-III surveys or precursors. Cluster samples are listed in
column 1. The references and the comparison surveys are listed in cols. 2, and
3, respectively. The number of clusters in the matched catalogue is reported
in col. 4. Differences between the reduced excess surface density, in
percents, are reported either considering the same weights for both the
reference and the comparison sample (col. 5), or using the respective survey
weights (col. 6).

Cluster sample | Survey ref. | Survey comp. | $N_{\rm cl}$ | weights ref | weights comp
---|---|---|---|---|---
CAMIRA-HSC-16A | HSC-SSP S16a | KiDS DR4 | $918$ |  $-9$$\,\pm\,$ | $10$ |  $-11$$\,\pm\,$ | $10$
AMICO-KiDS-DR3 | HSC-SSP S16a | KiDS DR4 | $1677$ |  $-6$$\,\pm\,$ | $10$ |  $-16$$\,\pm\,$ | $10$
AMICO-HSC-16A-XMM | HSC-SSP S16a | CFHTLenS | $2210$ |  $-3$$\,\pm\,$ | $17$ |  $-9$$\,\pm\,$ | $17$

Figure 12:  Average reduced excess surface density profiles of clusters
detected by CAMIRA in the fields covered by HSC-SSP S16a and KiDS DR4. Top
panel: average profiles. Bottom: differences.

Figure 13:  Same as Fig. 12 but for AMICO-KiDS-DR3 clusters.

Figure 14:  Average reduced excess surface density profiles of clusters
detected by AMICO in the XXL North field covered by HSC-SSP S16a and CFHTLenS.
Top panel: average profiles. Bottom: differences.

Shear calibration is regarded as one of the main systematic effects for
lensing and we can examine its effect by cross-checking shear estimates. We
compare measurements from overlapping survey areas. Direct comparison of
ellipticity estimates obtained with different resolution and/or in different
band passes might be misleading, as the spatial distribution of the light
emission may not be identical and different effective radial weight functions
may be used (Schrabback et al. 2018b). However, the estimated reduced
tangential cluster shear profile should be consistent when a matched shear
catalogue is used (Schrabback et al. 2018b).

In contrast with the main analysis in the paper, for this test we consider
cluster samples detected from the same data-sets used for shear measurements.
This maximises the relevant difference in signal. Source selection is
performed based on the photo-$z$ or colours of the reference survey. Any
difference in signal between the samples can be only due to differently
measured and calibrated galaxy shear.

We compare ellipticity measurements from a reference survey with a comparison
survey under two different schemes. We use either identical weights from the
reference survey, or we apply each respective survey weights consistently for
that survey. Results are summarised in Table 17. Agreement between HSC-SSP
S16a and KiDS DR4, see Figs. 12, 13, or CHFTLenS, see Fig. 14, is good, with
differences consistent with zero within the statistical uncertainty.


# Optimization of Sparse Convolution for 3D-Point Cloud on GPUs with CUDA

1st Chester Luo

Brion Technologies, an ASML Company
518000, Shenzhen, China
<EMAIL_ADDRESS>2nd Kevin Lai

School of Mechanical Engineering and Automation, Northeastern University
110000, Shenyang, China
<EMAIL_ADDRESS>
###### Abstract

The increasing utilisation of LiDAR and 3D sensors has made the analysis of 3D
point clouds indispensable in various applications, including object detection
and segmentation. In contrast to photos, point clouds possess distinctive
computing issues due to their sparse nature and absence of a regular grid.
Traditional Convolutional Neural Networks (CNNs) that are designed to perform
well on dense data are not suitable for this particular task. The utilisation
of sparse neural networks specifically engineered for the efficient handling
of sparse data is imperative in this context. Nevertheless, the utilisation of
neural networks on Graphics Processing Units (GPUs), particularly with the
CUDA framework, presents challenges due to the irregular nature of point cloud
data and the requirement for optimised memory access patterns. The present
study focuses on enhancing the efficiency of sparse convolution operators for
3D point clouds on GPUs through the utilisation of CUDA technology. This paper
presents a novel approach that combines the theoretical benefits of sparse
neural networks with efficient GPU-based implementations. By doing so, it
provides valuable insights and techniques for effectively using 3D point cloud
analysis, hence enhancing the capabilities of object detection and
segmentation in various fields.

###### Index Terms:

Sparse Convolution, CUDA, Point Cloud

##  I Introduction

In recent years, there has been a significant increase in the utilization of
deep learning methods, particularly convolutional neural networks (CNNs) [1],
which have emerged as the dominant approach in various domains that involve
structured grid data, such as picture analysis and processing. Nevertheless,
the exponential growth in the utilization of LiDAR and 3D sensors across many
domains has resulted in an increased need for the analysis of 3D point clouds.
The utilization of 3D point clouds is crucial in various applications,
including object recognition and segmentation, as they offer a spatial
depiction of things within a three-dimensional environment. In contrast to
photos, point clouds exhibit sparsity and lack a regular grid, hence posing
distinct processing and computational issues.

The examination of 3D point cloud data has initiated a novel period of
progress in the field of computer vision, enabling a wide range of
applications that use its abundant spatial information. One of the most
popular application of point cloud data is automated driving [2] and deep
learning can play role in various task on point cloud such as 3D shape
classification, 3D object detection and tracking, 3D point cloud segmentation
and 3D point cloud registration [3]. Qi et al. introduced a significant
contribution in the form of ”PointNet” [4], which marked a notable departure
from existing approaches. This study presented a deep learning framework that
was tailored to handle point sets, with a particular emphasis on addressing
challenges related to 3D classification and segmentation tasks. In addition,
Wang et al. introduced a real-time semantic segmentation method named
”PointSeg” [5] that focuses on utilizing 3D LiDAR point clouds for this
purpose. The aforementioned advancements highlight the extensive capabilities
of 3D point cloud analysis in diverse practical situations, ranging from
autonomous navigation to virtual reality.

Point clouds, distinct from traditional image data due to their sparsity,
present challenges in applying standard CNN architectures, which are better
suited for dense data. This sparsity, a result of the data collection
methodology and inherent information gaps in 3D space, necessitates the
adoption of sparse neural networks. These networks are tailored to efficiently
process sparse data while reducing computational demands. Sparse Convolution
Neural Network (SCNN) is regarded as solution to applied deep learning on
sparse data like point cloud. However, implementing them, particularly on
GPUs, is complex due to the irregular data structures and the need for
optimized memory access patterns. Additionally, deploying deep learning models
from servers to desktop or edge-computing devices remains challenging.
Recently, NVIDIA’s introduction of the TensorRT [6] inference framework has
shown promise in addressing these issues. TensorRT, which allows custom plugin
development, enhances the feasibility of deploying neural networks,
specifically those utilizing sparse convolution, on embedded systems.

This research investigates the challenges of implement sparse convolution
efficiently utilising GPUs with CUDA, to improve the speed of performing
inference on sparse convolution operators for 3D point clouds. The main goal
is to establish a correlation between the potential advantages of sparse
neural networks and their practical application in real-life scenarios, and we
will provide CUDA solutions for three categories of sparse convolution
operators. This paper is structured as follows: Section II will provide a
comparison between sparse convolution and traditional convolution to highlight
their distinct characteristics. Additionally, it will present the current
research on this topic. Section III presents a highly efficient method for
implementing sparse convolution on CUDA. Section IV outlines our discussions
and conclusions.

##  II Background and Motivation

###  II-A Structure and Represention for Point Cloud Data

Point cloud data, a fundamental representation in 3D data analysis, consists
of data points within a three-dimensional coordinate system, typically
generated by 3D scanning technologies. These points, representing the external
surfaces of objects, are characterized by their density, distribution, and
precision. Traditional point cloud representation primarily involves raw
coordinates, sometimes enriched with attributes like color, normal, and
reflectance. However, the unstructured nature of raw point clouds poses
challenges in computer vision and graphic because they are inherently complex
and irregular, consist of a massive number of points and contain noise and
outliers. Voxelization, the conversion of continuous point cloud data into a
discrete 3D grid, addresses these challenges. Each voxel in the grid
represents the presence or characteristics of points within its boundaries,
thereby regularizing the data for computational efficiency in storage,
compression, and processing. [7]

Figure 1: Voxelization of 3D points: (a) Original TLS points; and (b)
Voxelized point cloud. Source [7] The structured voxel grid offers a
simplified and computationally efficient approach for 3D data analysis,
suitable for applications in digital elevation modeling, urban planning, and
3D simulations.

When they need to be processed by neural network, for a point cloud and its
voxelized data, the sparse tensor might encapsulate the actual 3D coordinates
of points, alongside any associated attributes such as color, intensity, or
normal. In typical sparse convolution implementations, data are stored in
memory using specialized structures to efficiently represent and process the
sparsity. The indices structure holds the coordinates of non-zero points in a
sparse tensor, which corresponds to voxels or points in a point cloud. This
structure allows the algorithm to iterate only over the meaningful data
points, significantly reducing the computational load. For each batch in a
dataset, the indices are stored with their batch ID, allowing for parallel
processing of multiple data instances. The data structure for this work [8]
[9]

${0}$${x_{0}^{0}}$${y_{0}^{0}}$${z_{0}^{0}}$${0}$${x_{0}^{1}}$${y_{0}^{1}}$${z_{0}^{1}}$${\vdots}$${\vdots}$${\vdots}$${\vdots}$${1}$${x_{1}^{0}}$${y_{1}^{0}}$${z_{1}^{0}}$${1}$${x_{1}^{1}}$${y_{1}^{1}}$${z_{1}^{1}}$${\vdots}$${\vdots}$${\vdots}$${\vdots}$${M}$${x_{M}^{0}}$${y_{M}^{0}}$${z_{M}^{0}}$${M}$${x_{M}^{1}}$${y_{M}^{1}}$${z_{M}^{1}}$${\vdots}$${\vdots}$${\vdots}$${\vdots}$${M}$${x_{M}^{n}}$${y_{M}^{n}}$${z_{M}^{n}}$

$\vphantom{x_{1}^{1}}{D}^{1}_{1}$$\vphantom{x_{1}^{1}}{D}^{1}_{2}$$\vphantom{x_{1}^{1}}{D}^{1}_{3}$$\vphantom{x_{1}^{1}}\cdots$$\vphantom{x_{1}^{1}}{D}^{1}_{C}$$\vphantom{x_{1}^{1}}{D}^{2}_{1}$$\vphantom{x_{1}^{1}}{D}^{2}_{2}$$\vphantom{x_{1}^{1}}{D}^{2}_{3}$$\vphantom{x_{1}^{1}}{D}^{2}_{C}$$\vphantom{x_{1}^{1}}{D}^{3}_{1}$$\vphantom{x_{1}^{1}}{D}^{3}_{2}$$\vphantom{x_{1}^{1}}{D}^{3}_{3}$$\vphantom{x_{1}^{1}}{D}^{3}_{C}$$\vphantom{x_{1}^{1}}\vdots$$\vphantom{x_{1}^{1}}\ddots$$\vphantom{x_{1}^{1}}{D}^{N}_{1}$$\vphantom{x_{1}^{1}}{D}^{N}_{2}$$\vphantom{x_{1}^{1}}{D}^{N}_{3}$$\vphantom{x_{1}^{1}}{D}^{N}_{C}$$\vphantom{x_{1}^{1}}{D}^{1}_{1}$$\vphantom{x_{1}^{1}}{D}^{1}_{2}$$\vphantom{x_{1}^{1}}{D}^{1}_{3}$$\vphantom{x_{1}^{1}}\cdots$$\vphantom{x_{1}^{1}}{D}^{1}_{C}$$\vphantom{x_{1}^{1}}{D}^{2}_{1}$$\vphantom{x_{1}^{1}}{D}^{2}_{2}$$\vphantom{x_{1}^{1}}{D}^{2}_{3}$$\vphantom{x_{1}^{1}}{D}^{2}_{C}$$\vphantom{x_{1}^{1}}{D}^{3}_{1}$$\vphantom{x_{1}^{1}}{D}^{3}_{2}$$\vphantom{x_{1}^{1}}{D}^{3}_{3}$$\vphantom{x_{1}^{1}}{D}^{3}_{C}$$\vphantom{x_{1}^{1}}\vdots$$\vphantom{x_{1}^{1}}\ddots$$\vphantom{x_{1}^{1}}{D}^{N}_{1}$$\vphantom{x_{1}^{1}}{D}^{N}_{2}$$\vphantom{x_{1}^{1}}{D}^{N}_{3}$$\vphantom{x_{1}^{1}}{D}^{N}_{C}$$\vphantom{x_{1}^{1}}{D}^{1}_{1}$$\vphantom{x_{1}^{1}}{D}^{1}_{2}$$\vphantom{x_{1}^{1}}{D}^{1}_{3}$$\vphantom{x_{1}^{1}}\cdots$$\vphantom{x_{1}^{1}}{D}^{1}_{C}$$\vphantom{x_{1}^{1}}{D}^{2}_{1}$$\vphantom{x_{1}^{1}}{D}^{2}_{2}$$\vphantom{x_{1}^{1}}{D}^{2}_{3}$$\vphantom{x_{1}^{1}}{D}^{2}_{C}$$\vphantom{x_{1}^{1}}{D}^{3}_{1}$$\vphantom{x_{1}^{1}}{D}^{3}_{2}$$\vphantom{x_{1}^{1}}{D}^{3}_{3}$$\vphantom{x_{1}^{1}}{D}^{3}_{C}$$\vphantom{x_{1}^{1}}\vdots$$\vphantom{x_{1}^{1}}\ddots$$\vphantom{x_{1}^{1}}{D}^{N}_{1}$$\vphantom{x_{1}^{1}}{D}^{N}_{2}$$\vphantom{x_{1}^{1}}{D}^{N}_{3}$$\vphantom{x_{1}^{1}}{D}^{N}_{C}$BatchChannelsRows

Figure 2: Illustration of the data structures for sparse tensor for voxel or
point cloud. On the left is the indices structure, depicting the point
coordinates (m, x, y, z) for each batch ranging from 0 to M+1, where m is the
batch ID. On the right is the features structure, showcasing the organization
of feature values across channels, rows, and batches.

###  II-B Dense Convolution v.s. Sparse Convolution

In machine learning, ’sparse’ and ’dense’ matrices or tensors differ in their
non-zero element distribution. Dense matrices/tensors are predominantly non-
zero, implying rich information in almost every entry. In contrast, sparse
matrices/tensors have a substantial proportion of zeros (or a default value),
offering unique storage and computational challenges. Dense structures usually
require resources proportional to their size for storage and processing.
However, sparse structures benefit from specialized formats like Compressed
Sparse Row (CSR)[10] or ELLPACK Sparse Block (ESB)[11], optimizing non-zero
element storage and computation, thereby enhancing space efficiency.

3 1 1 1 0 2 0 2 0 1 1 1 0 0 1 2 0 1 1 1 14 3 0 1 1 0 03 0 1 1 0 2 40 1 1 0 4 0
51 1 0 0 5 0 02D dense matrix0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 10 0 0 0
1 0 00 0 1 0 0 0 00 0 1 0 0 0 01 1 0 0 0 0 02D sparse matrix Figure 3:
Illustrative comparison of a 2D dense matrix and a 2D sparse matrix. The dense
matrix predominantly features non-zero elements, whereas the sparse matrix
consists mainly of zero values with a few non-zero entries scattered
throughout. Such distinctions highlight the storage and computational
differences between the two matrix types.

In conventional convolutional procedures applied to dense data, a filter or
kernel traverses the input data, often an image, resulting in the generation
of an output feature map. The procedure is performing element-wise
multiplication between the kernel and a specific local region of the input,
and subsequently calculating the sum of the resulting products. In order to
enhance the efficiency of this computation, a commonly utilized technique
known as im2col is implemented. The im2col technique is employed to transform
localized areas of the input data into column vectors, effectively
transforming the convolution operation into a matrix multiplication, which is
especially friendly to GPU by using GEMM [12].

Figure 4: Illustration of the convolution process using the ‘im2col‘ approach,
source [13]. Starting with a set of kernels $K$ and an input $I$, the method
reshapes the input into overlapping patches. Correspondingly, kernels are
reshaped into rows to form the kernel-patch-matrix $K^{\prime}$. Matrix
multiplication between $K^{\prime}$ and the input-patch-matrix $I^{\prime}$
yields the output-patch-matrix $O^{\prime}$. The final output $O$ is derived
by reshaping $O^{\prime}$ to its intended dimensions.

Traditional convolution, though effective for dense data, struggles with
sparse datasets such as point clouds. Standard convolution assumes a dense,
consistent sampling, an assumption not applicable to the irregular, sparse
nature of point clouds. Consequently, these methods, including im2col, often
process extensive empty spaces in sparse data, resulting in increased
computational overhead and memory usage. To delve into the computational
complexity, consider a typical convolution operation that might require
$O(n^{2})$ operations for a dense $n\times n$ matrix. For sparse data with
only $k$ non-zero elements, an optimal approach would ideally require close to
$O(k)$ operations. However, with traditional convolution using techniques like
im2col, we still approach a complexity of $O(n^{2})$, even though a vast
majority of the processed elements are zeros. This means we’re expending
computational effort on a high proportion of elements that contribute no
meaningful information to the output. Such inefficiencies underline the
inadequacy of traditional convolution methods for sparse datasets and
emphasize the need for specialized techniques tailored to the intricacies of
sparse data.

Liu et al. introduced sparse convolutional neural networks (SCNN) in 2015
[14], an innovative approach leveraging a two-stage decomposition to reduce
inter-channel and intra-channel redundancies in kernel weights. The process
begins with an initial decomposition based on kernel weight reconstruction
error, followed by a fine-tuning stage incorporating a sparsity constraint.
This dual-phase approach effectively optimizes training error, kernel
sparsity, and the number of convolutional bases. Further delving into the SCNN
model, the authors elucidate that each sparse convolutional layer utilizes a
select few kernels followed by sparse matrix multiplication. While this
approach suggests inherent computational efficiency, the authors acknowledge
potential overheads in sparse matrix computations. To mitigate this, they
propose an optimized algorithm for sparse matrix multiplication, capitalizing
on the fixed nature of sparse kernels post-training to avoid indirect and
discontinuous memory access challenges. This CPU-focused strategy outperforms
standard sparse matrix libraries, offering significant acceleration over
traditional dense networks. This advancement is particularly beneficial for
deployment in typical clusters lacking specialized GPU resources.

In 2017, a new sparse convolution method named Submanifold[15] are introduced.
This type of convolution operates exclusively on the spatial locations where
input features are present, entirely bypassing the regions devoid of
information. SubM convolution leverages the spatial sparsity to reduce
computational requirements, ensuring that only relevant regions of the data
are processed. This efficiency makes it particularly suitable for point cloud
data where vast portions may lack relevant information. Based on this, Graham
et al. presented a novel approach to 3D semantic segmentation by employing
Submanifold Sparse Convolutional Networks [16] in 2018, show the effective
utilization of the distinctive characteristics of 3D data to achieve accurate.

0 1 1 1 0 0 00 0 1 1 1 0 00 0 0 1 1 1 00 0 0 1 1 0 00 0 1 1 0 0 00 1 1 0 0 0
01 1 0 0 0 0 0$\mathbf{I}$$*$1 0 1 0 1 0 1 0 1 $\mathbf{K}$$=$0 4 3 4 00 0 4 3
00 0 3 4 00 3 3 0 03 3 0 0 0$\mathbf{I*K}$$\times 1$$\times 0$$\times
1$$\times 0$$\times 1$$\times 0$$\times 1$$\times 0$$\times 1$ Figure 5:
Illustration of submanifold convolution on a 2D matrix. The input matrix $I$
is convolved with the kernel $K$. Only the central non-zero elements
(highlighted in orange) undergo computation. Regions with a central zero
(highlighted in red) are excluded from the convolution process.

###  II-C Related Work

In 2018, Elgammal et al. proposed a framework named uSCNN [17], an effective
sparse CNN inference engine that takes use of the sparsity of models and
activations, utilizing optimized sparse-sparse matrix convolution kernels with
compressed operands. This work optimized the sparse CNN on image. In the realm
for point cloud, Yan et al. introduces a general process to compute sparse
convolution includes building rules table and set input values for GEMM
process [18]. This work presents a comprehensive workflow detailing the entire
process of computation within the realm of sparse data from LiDAR on GPU.
Their work includes two stages: (1) Generate computation Rules (Algorithm 1).
(2) Compute values based on rules. The algorithm 1 desribes the whole process
and their code has been open-source on GitHub, names “SpConv” [19].

Algorithm 1 Rule Generation Process in [18]

0: Indice: Coordinates of active points, with dimensions of $n_{in}\times 3$;
Grid: A buffer with dimensions of $N\times D\times H\times W$; $i_{in}$: The
input index; $i_{out}$: The output index;

0: Rule

1: Rule[:] $\leftarrow-1$

2: $n_{in}\leftarrow$ number of input points

3: $n_{kernel}\leftarrow$ volume of convolution kernel

4: $n_{out}\leftarrow$ number of output points

5: for $i_{in}=0$ to $n_{in}$ do

6: $p_{in}\leftarrow$ Indice[$i_{in}$]

7: $p_{out}\leftarrow$ getOutputCoord($p_{in}$)

8: for p $\in p_{out}$ do

9: index $\leftarrow$ getSpatialIndex(p)

10: offset $\leftarrow$ getOffset(p, $p_{in}$)

11: Rule[offset, $i_{in}$, 0] $\leftarrow i_{in}$

12: Rule[offset, $i_{in}$, 1] $\leftarrow$ index

13: end for

14: end for

15: SpatialIndex $\leftarrow$ Rule[:, :, 1]

16: $n_{out}\leftarrow$ unique(SpatialIndex)

17: for $i_{out}=0$ to $n_{out}$ do

18: Grid[SpatialIndex[$i_{out}$]] $\leftarrow i_{out}$

19: end for

20: for $i_{in}=0$ to $n_{in}$ do

21: for $j=0$ to $n_{kernel}$ do

22: index $\leftarrow$ Rule[$j,i_{in},1$]

23: if index $>0$ then

24: Rule[$j,i_{in},1$] $\leftarrow$ Grid[index]

25: end if

26: end for

27: end for

In 2022, Tang et al. introduced a new point cloud inference engine named
“TorchSparse” [20], which has faster speed than ”SpConv”. This work introduces
innovative matrix multiplication grouping strategies that target efficiency in
computational tasks. The proposed fixed grouping strategy compensates for the
excess FLOPs with enhanced regularity, while the adaptive grouping seeks an
optimal balance point. Once the maps are devised, sparse convolution is
employed, further enhancing the computation efficiency. Significantly, the
utilization of matrix-vector multiplication on GPUs has been traditionally
low. Addressing this, the work integrates a gather-matmul-scatter computation
flow, ensuring a more streamlined and efficient process. Their work are also
open-source on Github [21]

Figure 6: TorchSparse concentrates on accelerating Sparse Convolution,
delineated into four key stages: mapping, gather, matrix multiplication, and
scatter-accumulate. Its objectives include augmenting the uniformity of sparse
workloads and reducing memory consumption. These aims are realized through the
use of adaptively batched matrix multiplication and the implementation of
quantized, vectorized, locality-aware scatter/gather operations. Furthermore,
kernel fusion in the mapping stage enhances memory efficiency.

However, both ”SpConv” and “TorchSparse” are specific designed for PyTorch
[22]. The design innovations in Spconv and TorchSparse, tailored to PyTorch’s
unique features and data types, have significantly enhanced this environment’s
capabilities. Yet, their PyTorch-centric approach, beneficial in some aspects,
also marks their fundamental limitations. These systems, designed specifically
for PyTorch, lack flexibility for use with other platforms or frameworks. This
issue is especially evident when deploying in varied computational
environments or integrating with non-PyTorch systems. Our research aims to
overcome these limitations by developing a version of sparse convolution that
maintains the computational efficiency of Spconv and TorchSparse while also
being compatible with a wider array of platforms and frameworks. This approach
broadens their applicability and increases deployment flexibility.

##  III Our approach

This section will present the methodology employed for designing operators to
maximize parallelism in CUDA, utilizing fundamental CUDA optimization
strategies. Parallelism in CUDA is attained through the simultaneous execution
of numerous threads, organized into blocks, across the multiple cores of a
GPU. In addition to the aforementioned considerations, it is imperative to
also prioritize the optimization of data load efficiency. Data locality refers
to the practice of reusing data that is already fetched into the faster memory
hierarchies, minimizing the need to access the slower global memory. In
contrast, shared memory plays a crucial role in facilitating inter-thread
communication and data exchange as it is accessible to all threads inside a
block[23]. By effectively utilizing both local and shared memory, it is
possible to enhance data access patterns, decrease memory latency, and lessen
dependence on the comparatively slower global memory.

###  III-A Input and Output Structure

The aforementioned publications discuss sparse tensor structures, which
consist of two main data tensors: one that represents the indices of the
points, and another that contains the related data for these points. This
structure is developer-friendly and efficient in PyTorch when constructing and
training networks. In order to incorporate this functionality into CUDA, it
becomes necessary to modify these tensors to conform to conventional CUDA data
types. The architecture of our operator utilises a collection of CUDA arrays
consisting of integers (int) to represent the indices, along with another set
of CUDA arrays for data (float). This approach effectively replaces the usual
data tensors.

A critical issue in this context is maintaining the correspondence between
input and output indices in typical sparse convolution (down-sampling) and
efficiently using this mapping in inverse convolution (up-sampling). Our
strategy emphasizes simplicity and efficiency. Instead of complicating the
mapping process, we use the indices as direct input variables. An additional
phase in the operator kernels then re-establishes their interconnection,
ensuring efficient operation and optimal performance.

In addition, to further improve the efficient of load input indices, we make a
little change in input format to achieve coalesced access. In coalesced
access, consecutive threads access consecutive memory locations, facilitating
a linear and contiguous memory access pattern. This approach not only
optimizes memory bandwidth utilization but also significantly bolsters kernel
performance. In our earlier design, indices were encompassed within a singular
CUDA int array, adhering to a layout defined as:

$\begin{array}[]{c}n_{0},x_{0},y_{0},z_{0},n_{1},x_{1},y_{1},z_{1},n_{2},x_{2},y_{2},z_{2},\dots\end{array}$

While this configuration is intuitively structured, it doesn’t cater to the
nuances of ensuring coalesced memory access during CUDA-based operations. To
address this, our enhanced design segregates the layout across four individual
int arrays. This new structure is delineated as:

$\displaystyle n_{0},n_{1},\dots\,|\,$ $\displaystyle x_{0},x_{1},\dots\,|$
$\displaystyle y_{0},y_{0},\dots\,|\,$ $\displaystyle z_{0},z_{0},\dots{}$

This configuration ensures an optimized memory access pattern, pivotal for
efficient CUDA processing.

OperatorInt:Indices (n,x,y,z)Int:Old Indices
(Inv)Float:FeaturesInt:Indices’(n’,x’,y’,z’)Float:Features’ Figure 7: Operator
design showcasing input-output structures. The operator takes in features
(Float:Features) and indices (Int:Indices (n,x,y,z)) along with old indices,
for inverse convolution (Int:Old Indices) and produces transformed features
(Float:Features’) and indices (Int:Indices’ (n’,x’,y’,z’)).

###  III-B Sampling Indices Mapping

As we get a list of points coordinates with no order, it is essential to
determine a rules to help us record these points and find a better
representation in computing process. Given a point’s position before
convolution, denoted as $(x,y,z)$, our objective is to establish its
relationship post-convolution. The following steps elucidate this process:

####  III-B1 Position Identification

Determine if the point is at the start of a convolution cell using the
divisibility condition:

$\displaystyle x\mod s$ $\displaystyle=0,y\mod s$ $\displaystyle=0,z\mod s$
$\displaystyle=0,$ (1)

where $s$ is the stride size.

####  III-B2 Start Cell Determination

If the point satisfies the above condition, it lies at the start of a cell.
Otherwise, ascertain the start cell’s coordinates as:

$\displaystyle startCell_{X}$
$\displaystyle=\left\lfloor\frac{x}{s}\right\rfloor\times s,$ (2)
$\displaystyle startCell_{Y}$
$\displaystyle=\left\lfloor\frac{y}{s}\right\rfloor\times s,$ $\displaystyle
startCell_{Z}$ $\displaystyle=\left\lfloor\frac{z}{s}\right\rfloor\times s.$

####  III-B3 Offsets Calculation

Calculate the point’s offset within its respective cell:

$\displaystyle offset_{X}$ $\displaystyle=x-startCellX,$ (3) $\displaystyle
offset_{Y}$ $\displaystyle=y-startCellY,$ $\displaystyle offset_{Z}$
$\displaystyle=z-startCellZ.$

These offsets represent the position of the point within the convolution
kernel.

####  III-B4 Position Mapping in the Convolved Space

The resulting position in the convolved space is determined as:

$\displaystyle o_{x}$ $\displaystyle=\frac{startCell_{X}}{s},o_{y}$
$\displaystyle=\frac{startCell_{Y}}{s},o_{z}$
$\displaystyle=\frac{startCell_{Z}}{s}.$ (4)

####  III-B5 Kernel Position Mapping

The point’s relative position within the kernel, considering its offset, is
represented as:

$kernel\\_offset=offset_{X}\times k^{2}+offset_{Y}\times k+offset_{Z},$ (5)

where $k$ denotes the kernel size.

###  III-C Submanifold Convolution

In our approach, the computation process for sparse convolution can be
concluded into three steps: create location table (LCT), create offset table
(OFT) and compute result. The Submanifold convolution and sparse convolution
has same steps in workflow, but there are little difference for inverse
convolution because its special input structure and compute rules. In this
part, we will introduce our methodology for Submanifold and others will be
introduced in next part. Our solutions for Submanifold conlution are mainly
follow the steps by works by Y. Yan, but we optimize the final computing
process and data structure in computing:

####  III-C1 Create Location Table (LCT)

In this step, we will scan the input indices to setup the relationship between
a point take partcipate in computation and its location in inputs.

To linearize the 3D coordinates $(x,y,z)$ into a single index, the
transformation is given by:

$\text{index}=x+y\times\text{max}_{x}+z\times\text{max}_{x}\times\text{max}_{y}$
(6)

The inverse transformation, which recovers the 3D coordinates from a given
index, is described as:

$\displaystyle z$
$\displaystyle=\left\lfloor\frac{\text{index}}{\text{max}_{x}\times\text{max}_{y}}\right\rfloor$
(7) index
$\displaystyle=\text{index}-z\times\text{max}_{x}\times\text{max}_{y}$ (8)
$\displaystyle y$
$\displaystyle=\left\lfloor\frac{\text{index}}{\text{max}_{x}}\right\rfloor$
(9) $\displaystyle x$ $\displaystyle=\text{index}\mod\text{max}_{x}$ (10)

In our work, The LCT can be stored in both CUDA Hash table [24] or large CUDA
array. For scenarios where the domain is relatively confined, allocating a
three-dimensional CUDA array of dimensions X, Y, and Z is advantageous as its
has better performance compared to CUDA hash table.

⬇

1__global__ void createLCTKernel(int *locationTable, const int *n_values,
const int *x_values, const int *y_values, const int *z_values, int N, int
max_x, int max_y) {

2 int idx = blockIdx.x * blockDim.x + threadIdx.x;

3 if (idx < N) {

4 int x = x_values[idx];

5 int y = y_values[idx];

6 int z = z_values[idx];

7 locationTable[getIndex(x, y, z, max_x, max_y)] = idx;

8 }

9}

Listing 1: CUDA kernel for location table creation

####  III-C2 Create Offset Table (OFT)

The offset table is a crucial element that guides convolution operations in
sparse data environments, representing the relationship between spatial
indices for efficient convolutions in sparse structures. The construction
processes for SubM and Spconv operators differ significantly. A notable
feature of the SubM convolution is its ability to preserve the congruity of
input and output feature map dimensions, facilitating uninterrupted
information flow. Identifying valid positions, central to the computation
kernel, is vital. The rule table’s creation depends on accurately defining
spatial neighborhoods, influenced by the kernel size. Each voxel or point
resides within a receptive field determined by the kernel dimensions, with the
kernel’s half-length defining its radial extent. Nested loops systematically
explore the kernel’s spatial range. At each point, the LCT is consulted to
check for voxel activity. Active voxels are recorded in the rule table with
their relative offsets, enhancing future data access efficiency.

$x_{0}$ $y_{0}$ $\vdots$ $\vdots$ $31$ $31$ $\vdots$ $\vdots$ $32$ $32$
$\vdots$ $\vdots$ $33$ $33$ $\vdots$ $\vdots$ Input Indices0 0 0 0 00 O 0 0 00
0 P 0 00 0 0 Q 0 0 0 0 0 0$o_{0}^{0}$ $\dots$ $o_{0}^{4}$ $\dots$
$o_{0}^{8}$$\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$-1 $\dots$ $I_{O}$
$\dots$ $I_{P}$$\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$$I_{O}$ $\dots$
$I_{P}$ $\dots$ $I_{Q}$$\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$$I_{P}$
$\dots$ $I_{Q}$ $\dots$ -1$\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$Result
OFT Figure 8: Visualization of the rule construction process for a submanifold
convolution. The ’Input Indices’ matrix represents the non-zero entries of the
input feature map. The process involves scanning these non-zero entries to
identify their valid neighboring non-zero entries. The ’Result OFT’ (Offset
Table) captures the spatial relations of these entries by setting the
corresponding locations within the table, thereby encoding the connectivity
patterns necessary for the convolution operation.

The convolution process for Spconv presents a more intricate challenge when
compared to SubM, particularly in the construction and management of the
Offset Table (OFT). This complexity arises due to the inherent nature of
Spconv, where the shapes of input and output feature maps aren’t identical,
introducing a variable dimensionality into the convolution operations. An
extra step are added in steps of building LOT, which will count the number of
down-sample points and get the shape of output tensor. In the output space,
one point may contributed by many points in origin input space, it raise a
predominant challenge that emerges during the Spconv process is managing and
mitigating CUDA race conditions. Specifically, the nature of parallelized GPU
operations, where multiple threads might attempt to access or modify shared
memory locations simultaneously, introduces potential conflicts

To handle this, we use The primary goal of the spconv_createLookupKernel is to
craft a mapping scheme – a ”Rule-Mapping Book” that correlates the input and
output spaces, incorporating the variances introduced by parameters like
kernel size and stride. Once this mapping is established, it’s imperative to
synchronize operations to ensure that the Rule-Mapping Book is fully populated
before subsequent operations, thereby ensuring data integrity and consistency.

This challenge is exemplified in the setUniquePoints kernel. The kernel’s
primary objective is to uniquely identify the spatial indices in the feature
map and populate the d_idxList and d_offsetTable accordingly. Given that
multiple threads might simultaneously operate on overlapping or neighboring
spatial regions, there’s a risk of threads overwriting each other’s data or
reading inconsistent data states, leading to inaccuracies in the resultant
OFT.

Algorithm 2 Normal Convolution Computation

0: Input: $inF,rT,w,outF,N,inC,outC,kS$ {F represents Features, C represents
Channels, rT is RulesTable, kS is kernel size}

0: Output: $outF$ after convolution

1: $kV\leftarrow kS^{3}$

2: for $i=0$ to $N-1$ do

3: for $o=0$ to $outC-1$ do

4: $s\leftarrow 0$

5: for $k=0$ to $kV-1$ do

6: $nIdx\leftarrow rT[i\cdot kV+k]$

7: if $nIdx\neq-1$ then

8: for $c=0$ to $inC-1$ do

9: $wIdx\leftarrow o\cdot kV\cdot inC+k\cdot inC+c$

10: $s\leftarrow s+inF[nIdx\cdot inC+c]\cdot w[wIdx]$

11: end for

12: end if

13: end for

14: $outF[i\cdot outC+o]\leftarrow s$

15: end for

16: end for

####  III-C3 Compute Result

This step will execute computation process for output using OFT and weights.
The original method, as described in Algorithm 2, adopts a nested loop
approach to compute the convolution values. While conceptually
straightforward, this technique exhibits inefficiencies when implemented on a
GPU architecture. One of the main limitations is the frequent access to the
global memory, especially for retrieving weight values. Since global memory
access is notably slower than shared or local memory access on CUDA devices,
this becomes a performance bottleneck. The nested nature of the loops
amplifies this inefficiency, leading to suboptimal performance, especially for
larger volumes of data.

$o_{1}^{1}$ $o_{1}^{2}$ $o_{1}^{3}$ $\cdots$ $o_{1}^{c}$ $o_{2}^{1}$
$o_{2}^{2}$ $o_{2}^{3}$ $\cdots$ $o_{2}^{c}$ $o_{3}^{1}$ $o_{3}^{2}$
$o_{3}^{3}$ $\cdots$ $o_{3}^{c}$ $o_{4}^{1}$ $o_{4}^{2}$ $o_{4}^{3}$ $\cdots$
$o_{4}^{c}$ $o_{5}^{1}$ $o_{5}^{2}$ $o_{5}^{3}$ $\cdots$ $o_{5}^{c}$ $\vdots$
$\vdots$ $\vdots$ $\vdots$ $\vdots$ $o_{n}^{1}$ $o_{n}^{2}$ $o_{n}^{3}$
$\cdots$ $o_{n}^{c}$ Offset Table
(Features)$*$$w_{1}$$w_{2}$$w_{3}$$\vdots$$w_{n}$Weights$=$$r_{1}$$r_{2}$$r_{3}$$r_{4}$$r_{5}$$\vdots$$r_{n}$Results

$r_{1}^{1}$ $r_{1}^{2}$ $r_{1}^{3}$ $\cdots$ $r_{1}^{c}$ $r_{2}^{1}$
$r_{2}^{2}$ $r_{2}^{3}$ $\cdots$ $r_{2}^{c}$ $r_{3}^{1}$ $r_{3}^{2}$
$r_{3}^{3}$ $\cdots$ $r_{3}^{c}$ $r_{4}^{1}$ $r_{4}^{2}$ $r_{4}^{3}$ $\cdots$
$r_{4}^{c}$ $r_{5}^{1}$ $r_{5}^{2}$ $r_{5}^{3}$ $\cdots$ $r_{5}^{c}$ $\vdots$
$\vdots$ $\vdots$ $\vdots$ $\vdots$ $r_{n}^{1}$ $r_{n}^{2}$ $r_{n}^{3}$
$\cdots$ $r_{n}^{c}$ $grid.y$$grid.x$

Figure 9: The upper portion illustrates the computation of convolution from a
given set of weights and a offset table. The offset table, representing
features, directs the mapping of weights to corresponding feature values by
indexed rules, facilitating the convolution process. The lower portion depicts
the CUDA grid configuration and the allocation of weights in shared memory.
Each output channel is assigned to a distinct CUDA grid dimension, with its
weights stored in shared memory to expedite the convolution operation.

The optimized method offers a more efficient approach. By assigning the output
channel, $oc$, directly to the block’s y-dimension, the nested loop for output
channels is eliminated. Furthermore, the most significant enhancement is the
introduction of shared memory, which caches the ‘weights‘ for faster access.
This approach dramatically reduces the frequency of slow global memory
accesses. Loading the ‘weights‘ into shared memory ensures that all threads
within a block can access the values with reduced latency, speeding up the
convolution computation considerably.

⬇

1__global__ void computeConv(const float *input,

2 const int *rules, const float *w, float *output,

3 int num, int in_ch, int out_ch, int k_size) {

4 int idx = blockIdx.x * blockDim.x + threadIdx.x;

5 int oc = blockIdx.y;

6 int k_vol = k_size * k_size * k_size;

7 extern __shared__ float sW[];

8 for (int k = 0; k < k_vol; k++) {

9 for (int ic = 0; ic < in_ch; ic++) {

10 int w_idx = oc * k_vol * in_ch + k * in_ch + ic;

11 sW[k * in_ch + ic] = w[w_idx];

12 }

13 }

14 __syncthreads();

15

16 if (idx < num) {

17 float sum = 0;

18 for (int k = 0; k < k_vol; k++) {

19 int nIdx = rules[idx * k_vol + k];

20 if (nIdx != -1) {

21 for (int ic = 0; ic < in_ch; ic++) {

22 sum += input[nIdx * in_ch + ic] * sW[k * in_ch + ic];

23 }

24 }

25 }

26 output[idx * out_ch + oc] = sum;

27 }

28}

###  III-D Solution for Normal Sparse Convolution

As there are might many position map to a same point after down sample
process, there are two crucial aspects we need to consider in design sparse
convolution operator: (1) How to get the amount of elements after down sample,
which determines the output size and shape of offset table. (2) Build suitable
rules table for them. To handle this, we used a counter to record the output
size by combination of CUDA atomic opeation and a temp table to record a
position after downsample exists, and its location in output list. And for
each location in origin input, we used a new-design two-stage offset table.
The first table has only two column. The first column records index in origin
input second column records its offset in kernel. In the second stage,
according to above table, then we can set offset table with same schema in
above part used for convolution computing.

###  III-E Design for Inverse Convolution

The process of deconvolution, which is the mathematical inverse of
convolution, plays a crucial role in the technique of up-sampling. Up-sampling
is a method used to enhance the data rate or resolution. By employing the
technique of deconvolution, it becomes possible to recreate a dataset with
higher resolution by reversing the convolutional procedure that was applied to
its downsampled version. Fig shows process of inverse convolution.

Figure 10: Schematic diagram of convolution and inverse convolution process,
source: [25]. The inverse convolution essentially redistributes the values
from O across a larger matrix using the filter D, aiming to approximate the
original matrix X’s structure and values in Y.

In the proposed methodology, we introduce an innovative algorithm for inverse
convolution. The algorithm begins by considering two sets of indices: the
first (indices1) represents the data prior to the downsample process, and the
second (indices2) encapsulates the data after downsampling. This dual-index
approach forms the foundation for our offset table, crucial for the inverse
process. The core of our methodology lies in reversing the downsampling
coordinates computation. By meticulously scanning through indices1, we locate
the corresponding coordinates in indices2. This enables us to accurately
determine the offset in the original convolution kernel, a critical step for
setting the input feature values correctly and performing the convolution
operation effectively.

Our implementation is structured in three distinct stages: (1) Generate
location table for indices2. (2) Utilize the location table to correctly
assign input feature values to their appropriate locations in the temporary
feature array. (3) Perform the inverse convolution operation. It calculates
the output features by aggregating the weighted sum of input features and
weights, as defined by the convolution kernel.

##  IV Conclusion

We mainly contribute two improvement in the implement of sparse convolution:
(1) Optimize convolution computation process utlizing the advantage of CUDA.
(2) A easy and efficient approch to implement inverse sparse convolution. This
work effectively presents an innovative implementation of sparse convolution
operators using CUDA, emphasizing maximized parallelism and efficient data
load optimization. It introduces a novel approach to handling tensor
structures within CUDA’s framework, aligning with PyTorch’s usability while
leveraging CUDA’s performance capabilities. The methodology, encompassing the
creation of location table and offset table followed by computation, optimizes
sparse convolution processes, particularly addressing inefficiencies in
traditional nested loop approaches. By incorporating shared memory for weight
caching and optimizing thread and block utilization, the thesis demonstrates
significant advancements in the efficiency and speed of sparse convolution
operations on GPU architectures, offering substantial contributions to the
field of GPU-accelerated computing.

## References

  * [1] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _Nature_ , vol. 521, pp. 436–44, 05 2015.
  * [2] D. Fernandes, A. Silva, R. Névoa, C. Simões, D. Gonzalez, M. Guevara, P. Novais, J. Monteiro, and P. Melo-Pinto, “Point-cloud based 3d object detection and classification methods for self-driving applications: A survey and taxonomy,” _Information Fusion_ , vol. 68, pp. 161–191, 2021.
  * [3] Y. Guo, H. Wang, Q. Hu, H. Liu, L. Liu, and M. Bennamoun, “Deep learning for 3d point clouds: A survey,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 43, no. 12, pp. 4338–4364, 2021.
  * [4] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” 2017.
  * [5] Y. Wang, T. Shi, P. Yun, L. Tai, and M. Liu, “Pointseg: Real-time semantic segmentation based on 3d lidar point cloud,” _arXiv preprint arXiv:1807.06288_ , 2018.
  * [6] H. Vanholder, “Efficient inference with tensorrt,” in _GPU Technology Conference_ , vol. 1, no. 2, 2016.
  * [7] Y. Xu, X. Tong, and U. Stilla, “Voxel-based representation of 3d point clouds: Methods, applications, and its potential use in the construction industry,” _Automation in Construction_ , vol. 126, p. 103675, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0926580521001266
  * [8] Z. Liu, H. Tang, Y. Lin, and S. Han, “Point-voxel cnn for efficient 3d deep learning,” in _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds., vol. 32. Curran Associates, Inc., 2019.
  * [9] J. Deng, S. Shi, P. Li, W. Zhou, Y. Zhang, and H. Li, “Voxel r-cnn: Towards high performance voxel-based 3d object detection,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, pp. 1201–1209, 05 2021.
  * [10] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. Van der Vorst, _Templates for the solution of linear systems: building blocks for iterative methods_. SIAM, 1994.
  * [11] X. Liu, M. Smelyanskiy, E. Chow, and P. Dubey, “Efficient sparse matrix-vector multiplication on x86-based many-core processors,” in _Proceedings of the 27th International ACM Conference on International Conference on Supercomputing_ , ser. ICS ’13. New York, NY, USA: Association for Computing Machinery, 2013, p. 273–282. [Online]. Available: https://doi.org/10.1145/2464996.2465013
  * [12] “Matrix Multiplication Background User Guide - NVIDIA Docs — docs.nvidia.com,” https://docs.nvidia.com/deeplearning/performance/dl-performance-matrix-multiplication/index.html, [Accessed 01-11-2023].
  * [13] A. Vasudevan, A. Anderson, and D. Gregg, “Parallel multi channel convolution using general matrix multiplication,” in _2017 IEEE 28th international conference on application-specific systems, architectures and processors (ASAP)_. IEEE, 2017, pp. 19–24.
  * [14] B. Liu, M. Wang, H. Foroosh, M. Tappen, and M. Pensky, “Sparse convolutional neural networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 806–814.
  * [15] B. Graham and L. van der Maaten, “Submanifold sparse convolutional networks,” _CoRR_ , vol. abs/1706.01307, 2017. [Online]. Available: http://arxiv.org/abs/1706.01307
  * [16] B. Graham, M. Engelcke, and L. Van Der Maaten, “3d semantic segmentation with submanifold sparse convolutional networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 9224–9232.
  * [17] M. A. Elgammal, O. M. Awad, I. E. Vivancos, A. Moshovos, and V. Betz, “cuscnn: an efficient cuda implementation of sparse cnns,” in _Proceedings of the 13th International Symposium on Highly Efficient Accelerators and Reconfigurable Technologies_ , 2023, pp. 107–113.
  * [18] Y. Yan, Y. Mao, and B. Li, “Second: Sparsely embedded convolutional detection,” _Sensors_ , vol. 18, no. 10, 2018. [Online]. Available: https://www.mdpi.com/1424-8220/18/10/3337
  * [19] “GitHub - traveller59/spconv: Spatial Sparse Convolution Library — github.com,” https://github.com/traveller59/spconv, [Accessed 30-10-2023].
  * [20] H. Tang, Z. Liu, X. Li, Y. Lin, and S. Han, “Torchsparse: Efficient point cloud inference engine,” _Proceedings of Machine Learning and Systems_ , vol. 4, pp. 302–315, 2022.
  * [21] “GitHub - mit-han-lab/torchsparse: [MLSys’22] TorchSparse: Efficient Point Cloud Inference Engine — github.com,” https://github.com/mit-han-lab/torchsparse, [Accessed 01-11-2023].
  * [22] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in _Advances in Neural Information Processing Systems 32_. Curran Associates, Inc., 2019, pp. 8024–8035. [Online]. Available: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * [23] M. Garland, S. Le Grand, J. Nickolls, J. Anderson, J. Hardwick, S. Morton, E. Phillips, Y. Zhang, and V. Volkov, “Parallel computing experiences with cuda,” _IEEE micro_ , vol. 28, no. 4, pp. 13–27, 2008.
  * [24] D. A. Alcantara, A. Sharf, F. Abbasinejad, S. Sengupta, M. Mitzenmacher, J. D. Owens, and N. Amenta, “Real-time parallel hashing on the gpu,” in _ACM SIGGRAPH asia 2009 papers_ , 2009, pp. 1–9.
  * [25] S.-F. Wang, W.-K. Yu, and Y.-X. Li, “Multi-wavelet residual dense convolutional neural network for image denoising,” _IEEE Access_ , vol. 8, 11 2020.