TL;DR
This appears unlikely to work at scale, at least with current parameter choices. There's a classical post-processing step which is worse than the quadratic sieve in complexity unless the quantum part of the algorithm starts returning answers that are significantly better than those produced in initial experiments.
I'll try and walk through what's going on step-by-step.
Factoring 101: difference of squares
The last step of the method is well-understood and is the basis of almost all modern factoring methods. To factor a number $N$ we try and find two congruent squares modulo $N$:
$$a^2\equiv b^2\pmod N\iff a^2-b^2\equiv 0\pmod N.$$
Given two such integers $a$, $b$, we can use the high school identity $a^2-b^2=(a-b)(a+b)$ to deduce that
$$(a+b)(a-b)\equiv 0\pmod N$$
and so every prime divisor of $N$ divides either $a+b$ or $a-b$. If the prime divisors divide different brackets (which they will a positive proportion of the time) we can take $\mathrm{GCD}(a+b,N)$ and $\mathrm{GCD}(a-b,N)$ using Euclid's algorithm and recover factors.
Index Calculus: Making squares from smooth numbers
Modern factoring methods approach the challenge of finding congruent squares by making squares from products of smooth numbers. These are integers all of whose prime factors lie below some bound $B$. We can write smooth numbers in index notation which is a $\pi(B)+1$ long vector which records the exponents of the prime factorisation of the number. For example if $B=10$ then smooth numbers are those whose only prime divisors are 2, 3, 5 or 7 (and also possibly the unit -1), and I can write -700 in index notation thus:
$$-700=-1^1\times 2^2\times 3^0\times 5^2\times 7=(1,2,0,2,1).$$
Note that if I multiply two smooth numbers, I get a smooth number whose index representation is the sum of the two index representations. Also note that if all of the components of an index representation are even, this is precisely equivalent to a smooth number being a perfect square. This allows us to find products of smooth numbers that are squares using linear algebra on index representations: if we find more than $\pi(B)+1$ $B$-smooth numbers then there is a linear combination of their index representations all of whose entries are 0 mod 2 and we can find this linear combination e.g. with Gaussian elimination. Dixon's random squares method is a simple example of a factorisation method using this idea. The general number field sieve (the best known classical factoring algorithm) essentially uses a slightly more complex version of this idea.
There's a trade-off in the choice of $B$: choose $B$ too small and smooth number become rarer and very hard to find; choose $B$ too large and we have to find too many smooth numbers.
Schnorr: Finding smooth numbers using short vector methods
The previous two sections are well-established and quite well-understood, we now head into choppier waters. Established factoring algorithms have found the smooth numbers necessary to make squares by taking a large set of small candidates (small numbers are more likely to be smooth) and testing each of them for smoothness using souped-up trial division methods (sieving). Schnorr's lattice-factoring method aims to build a machine that spits out smooth numbers thereby massively saving on the search time. This also could allow us to use much smaller $B$ values.
Schnorr's idea is to find triples of smooth numbers $u$, $v$ such that $u-vN$ is very small indeed (and hence very likely to be smooth). This gives us
$$u\equiv u-vN\pmod N$$
and index calculus can then find a subset of the triples such that the product of both the $u$ and the $u-vN$ are both squares and we're home. The new paper gives examples of the triples from their smaller experiments on pages 22 and 24 of their paper.
Schnorr's idea is to write $u=p_0^{e_0}\cdots p_k^{e_k}$ and $v=p_0^{f_0}\cdots p_k^{f_k}$ and then take real logarithms to write
$$\log u=e_0\log p_0+\cdots +e_k\log p_k$$
$$\log vN=f_0\log p_0+\cdots +f_k\log p_k+\log N.$$
The idea is that if $|u-vN|$ is small, then $u\approx vN$ and so $\log u\approx\log(vN)$ and so $\log u-\log(vN)$ is small. Moreover $\log u-\log(vN)$ is an integer linear combination of the known real numbers $\{\log p_i:p_i<B\}$ and $\log N$ and spitting out smooth values $u$ and $v$ should be achievable by spitting out small linear combinations of logarithms. Finding small linear combinations is a problem that can be expressed in terms of a close vector problem. If we consider the lattice of integer linear combinations of rows of the matrix
$$\begin{pmatrix}1 & 0 & 0 &\cdots & 0 & 0 & \log p_0\\
0 & 1 & 0 &\cdots & 0 & 0 & \log p_1\\
0 & 0 & 1 &\cdots & 0 & 0 & \log p_2\\
\vdots & & & \ddots & & \vdots &\vdots\\
0 & 0 & 0 &\cdots & 1 & 0 & \log p_{k-1}\\
0 & 0 & 0 & \cdots & 0 & 1 & \log p_k\\
\end{pmatrix}$$
and we find a linear combination close to the vector $(0\ 0\ 0\ldots 0\ 0\ \log N)$, then we have found a pair of reasonably-sized smooth numbers $u$ and $v$ such that $\log u-\log(vN)$ is small. (Schnorr's method gives different weights to the columns to produce many solutions and also to try to guarantee that $u-vN$ is really small; you can think of this as reducing the above lattice with respect to a weighted Euclidean metric).
The idea is compelling, but there are issues*. Let us consider how close our close vector needs to be. As we will be naively testing $u-vN$ for smoothness, we should aim to find $u$ and $v$ values where $u-vN$ is significantly smaller than existing methods such as the quadratic sieve (where we test number of size $N^{1/2}$) or number field sieve (where we test numbers of size $\exp(c(\log N)^{2/3}(\log\log N)^{1/2})<N^{o(1)}$. If the final component of the difference vector in our close vector solution is $\epsilon$ then we have a solution where
$$\log u=\log(vN)+\epsilon\iff u=vN\exp(\epsilon)$$
so that
$$u-vN=vN(\exp(\epsilon-1)>\epsilon vN.$$
For this to be comparable to the quadratic sieve we will need $\epsilon$ to be less than $N^{-1/2}$ and to be competitive with the number field sieve it needs to be less than $N^{-1+o(1)}$. These are incredibly small which not only puts a strain on lattice algorithm, but also on the accuracy of our inputs: each of the logarithms of small primes will need to be expressed to 100s of decimal places of accuracy else the approximation will swamp $\epsilon$. This is something to watch closely as problems scale up.
QAOA: Finding short vectors using quantum computers
The novelty within the new paper seems to be a use of QAOA (quantum approximate optimisation algorithm/ansatz) to improve solutions to close vector problems. I've not yet read anything in the paper that makes their idea unique to the Schnorr lattice. The idea is that given a close-ish vector $\mathbf b=(b_0\ b_1\ldots b_k)$ to a target vector $\mathbf t=(t_0\ t_1\ldots t_k)$ they use QAOA to search the region
$$\{(\mathbf v=(v_0\ v_1\ldots v_k):|v_i-b_i|\le 1\}$$
(i.e. the hypercube of side length 2 centred on our close-ish vector) for an approximate minimum to the function**
$$||\mathbf v-\mathbf t||_2.$$
The minimum of this function would be taken at the closest vector to $\mathbf t$ in the hypercube. This is a search space of size $3^{k+1}$ and so brute force exhaustive search is not a sensible approach for large dimensions. The factoring application will rely on our really close vector lying within the hypercube, which will depend on the starting point. The paper says to derive starting points from Babai's rounding algorithm (i.e. rounding coordinate-wise with respect to the basis) applied to an LLL reduced basis. This will not provide spectacularly good starting points for higher dimensions, but it may be that the process can be iterated to "hill-climb" to iteratively better solutions. As written however, the paper relies on classical computation having found a close-ish vector that has a really close vector within its hypercube neighbourhood and this feels like a big assumption as we scale.
How good is QAOA at finding approximate minima to multivariate functions in hypercubes? Answers to this question are arcane. The behaviour of the function in the hypercube is modelled as a Hamiltonian, this Hamiltonian is modelled in a circuit and then the energy of the Hamiltonian approximately minimised. Behaviour varies according to how well-conditioned the Hamiltonian is and leaving questions of eigenvalues and Pauli states to quantum.stackexchange.com, I'll say that it's certainly likely to find something closer than the initial Babai starting point, though how much closer I can't say.
Issues with the experimental data
My increasing scepticism derives from looking at the quoted data. Schnorr's method spits out smooth $u$ and $v$ values for which $u-vN$ is supposed to be really small, so that it is automatically smooth. For larger values of $u-vN$ we have to do naive smoothness tests much as we do for existing factoring algorithms and if the $u-vN$ values are around $\sqrt N$ in size then the whole method is unlikely to beat the quadratic sieve, let alone the number field sieve. Looking at p. 11 (Section III.A. of the supplementary) they quote a Lemma of Schnorr*** that says $u-vN$ values of size around $O(p_k)$ should be produced. These would certainly be suitable and very likely to be smooth. The authors conclude Section III by relying on this estimate. In all of their experiments though, $u-vN$ comes up much larger than $p_k$.
In their "3 qubit" experiment to factor 1961 their (scaled) target is $(0\ 0\ 0\ 240)$ and an example Babai starting point of $(0\ 4 \ 4\ 242)$ which corresponds to $u=2025=3^45^2$ $v=1$ and a $u-vN$ value of 64 (which although smooth is considerably bigger than $p_k=5$). For this size of problem it is possible to exhaust the hypercube and note that there are no closer vectors. On page 22 they list all of the factorisations of suitable $u-vN$ that they found. These range in size from 17 to 5453. All of them are bigger than $p_k$, all bar three are bigger than $\sqrt N$ and some are bigger than the number that we wish to factor.
Likewise for their "5 qubit" case (factoring 48567227) on page 24 the value of $u-vN$ varies between 256476 and 222890309. All of them have a prime factor bigger than 11 and all of them are considerably bigger than $\sqrt N$.
In the "10 qubit" case (factoring 261980999226229) on page 19 they quote example $u-vN$ values between 10 and 11 decimal digits all bigger than $N^{0.7}$.
They rescue these examples by declaring a different smoothness bound for the $u-vN$ and naively testing them for smoothness (See p. 16 Section IV.D), Such an approach will be worse in classical post-processing than the quadratic sieve unless the $u-vN$ value can be made consistently less than $N^{1/2}$. Early experimental data does not indicate that this is the case.
(*) - The transcendental number theorist in me is queasy because there are deep results saying that linear forms in logarithms cannot be too small and that finding instances of smooth plus smooth equals smooth are tricky, but I'll try to stick to more concrete worries.
(**) - For the factoring application, a better approach would be simply to minimise the value of the final component which is all we really care about.
(***) - I haven't seen a satisfactory proof of this lemma.
