@lcz Very nice. Looking through the paper, it seems that rather than take the pseudoinverse of L, you can take the regular inverse of L+(1/n)*J, where J is the matrix of 1's. (The wikipedia page seems to have missed this point.) Since you are not making use of GraphTheory's LaplacianMatrix, and making your own weighted version, you may as well generate the above matrix directly (called LJ in the code). Its diagonals are chosen so the row/column sums are 1. I'm assuming the regular inverse is more efficiently calculated than the Moore-Penrose pseudoinverse, though I didn't check. This leads to the code below (works for weighted and unweighted graphs).
if IsDirected(G) then error "Must be an undirected graph" end if;
if IsWeighted(G) then W:=WeightMatrix(G) else W:=WeightMatrix(MakeWeighted(G)) end if;
M1:=Matrix(n,n, (i,j)-> if W[i,j]<>0 then 1/n-1/W[i,j] else 1/n end if,'shape'='symmetric');