In recent years, substantial progress has been made on graph convolutional networks (GCN). In this paper, for the first time, we theoretically analyze the connections between GCN and matrix factorization (MF), and unify GCN as matrix factorization with co-training and unitization. Moreover, under the guidance of this theoretical analysis, we propose an alternative model to GCN named Cotraining and Unitized Matrix Factorization (CUMF). The correctness of our analysis is verified by thorough experiments. The experimental results show that CUMF achieves similar or superior performances compared to GCN. In addition, CUMF inherits the benefits of MF-based methods to naturally support constructing mini-batches, and is more friendly to distributed computing comparing with GCN. The distributed CUMF on semi-supervised node classification significantly outperforms distributed GCN methods. Thus, CUMF greatly benefits large scale and complex real-world applications.
The work poses an interesting question: Are GCNs (and GNNs) just special types of matrix factorization methods? Unfortunately, the short answer is **no**, which goes against what the authors say. 
Until recently I thought like the authors, but the concurrent work [1] (On the Equivalence between Node Embeddings and Structural Graph Representations) https://openreview.net/forum?id=SJxzFySKwH changed my mind. 
The work of (Li et al., 2018) shows that nearby nodes tend to get similar representations. There is mounting experimental evidence of that being the case in real-world graphs (e.g., https://arxiv.org/abs/1908.08572). But [1] shows that GCNs and GNNs are fundamentally different from matrix factorization methods, regardless of the loss function used to learn the embeddings. Consider Figure 1 in [1], and it is easy to see that matrix factorization will give different embeddings to the Lynx and Orca nodes, while GCNs and GNNs must give the same embedding. Even if we connect the graphs through the Spruce and the Zooplankton nodes, their conclusion would not change. Matrix factorization (as broadly understood) will give embeddings that can even be used to cluster nodes. The eigenvectors of the symmetric Laplacian encode the diffusion of a type of random walk and nodes that are far away in the graph must have different embeddings (because through the diffusion operator, they are far away).
In GCNs, the convergence of the embeddings is better explained by the mixing of a random walk (Theorem 1 of (Xu et al., 2018)), which, in the special case of a GCN, converges to 1/sqrt(degree of node), as shown by (Li et al., 2018) in their Theorem 1 for the symmetric Laplacian. This is unrelated to what we get in matrix factorization as explained earlier. 
What is wrong with the math: Equation (11) is equated with matrix factorization, but note that it does not account for nonedges, while matrix factorization accounts for nonedges. This issue is more clear in Equation (14). The problem happens when the paper jumps from Equation (14), which is correct but not MF, to Equation (15) which is MF but unrelated to Equation (14). The argument is that “negative edges sampling is used, for better convergence”… sorry, not for better convergence, it completely changes the optimization objective. Hence, GCNs are not matrix factorization methods.
I think the paper is a valiant effort, but unfortunately the core premise is incorrect. The jump from Equation (14) to equation (15) cannot be justified, and I believe showcases a fundamental flaw the argument. I do not see a way to fix the paper. I vote to reject it. 
[1] On the Equivalence between Node Embeddings and Structural Graph Representations, https://openreview.net/forum?id=SJxzFySKwH 
Xu, K., Li, C., Tian, Y., Sonobe, T., Kawarabayashi, K.I. and Jegelka, S., 2018. Representation learning on graphs with jumping knowledge networks. ICML 2018.
Li, Qimai, Zhichao Han, and Xiao-Ming Wu. Deeper insights into graph convolutional networks for semi-supervised learning. AAAI, 2018.
--------------
Read rebuttal. Will keep my original assessment.