{'fraction_non_alphanumeric': 0.04010273624249822, 'fraction_numerical': 0.023268233140067134, 'mean_word_length': 4.010830784913353, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 39, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing low-rank approximations of large-scale matrices with the Tensor Network randomized SVD

We propose a new algorithm for the computation of a singular value decomposition (SVD) low-rank approximation of a matrix in the Matrix Product Operator (MPO) format, also called the Tensor Train Matrix format. Our tensor network randomized SVD (TNrSVD) algorithm is an MPO implementation of the randomized SVD algorithm that is able to compute dominant singular values and their corresponding singular vectors. In contrast to the state-of-the-art tensor-based alternating least squares SVD (ALS-SVD) and modified alternating least squares SVD (MALS-SVD) matrix approximation methods, TNrSVD can be up to 17 times faster while achieving the same accuracy. In addition, our TNrSVD algorithm also produces accurate approximations in particular cases where both ALS-SVD and MALS-SVD fail to converge. We also propose a new algorithm for the fast conversion of a sparse matrix into its corresponding MPO form, which is up to 509 times faster than the standard Tensor Train SVD (TT-SVD) method while achieving machine precision accuracy. The efficiency and accuracy of both algorithms are demonstrated in numerical experiments.

1. Introduction. When Beltrami established the existence of the singular value decomposition (SVD) in 1873 [2,20], he probably had not foreseen that this matrix factorization would become a crucial tool in scientific computing and data analysis [3,9,10,13]. Among the many applications of the SVD are the determination of the numerical rank and condition number of a matrix, the computation of low-rank approximations and pseudoinverses, and solving linear systems. These applications have found widespread usage in many fields of science and engineering [10,20]. Matrices with low numerical ranks appear in a wide variety of scientific applications [15]. For these matrices, finding a low-rank approximation allows them to be stored inexpensively without much loss of accuracy. It is not uncommon for matrices in data analysis to be very large and classical methods to compute the SVD [11,12,13] can be ill-suited to handle large matrices. One proposed solution to compute a low-rank approximation of large data matrices is to use randomized algorithms [15]. An attractive feature of these algorithms is that they require only a constant number of passes over the data. It is even possible to find a matrix approximation with a single pass over the data [14]. This enables the efficient computation of a low-rank approximation of dense matrices that cannot be stored completely in fast memory [30].
Another way to handle large matrices is to use a different data storage representation. A Matrix Product Operator (MPO), also called Tensor Train Matrix [24], is a particular tensor network representation of a matrix that originates from the simulation of one-dimensional quantum-lattice systems [29]. This representation transforms the storage complexity of an n d × n d matrix into O(dn 2 r 2 ), where r is the maximal MPO-rank, effectively transforming the exponential dependence on d into a linear one. The main idea of this representation is to replace the storage of a particular matrix element by a product of small matrices. An efficient representation is then found when only a few small matrices are required. A Matrix Product State (MPS), also called a Tensor Train [22], is a similar tensor network representation of a vector. These tensor network representations have gained more interest over the past decade, together with their application to various problems [1,4,5,6,19,21,25,26]. In particular, finding an SVD low-rank approximation of a matrix in the MPO representation is addressed in [18] with the development of the ALS-SVD and MALS-SVD methods. These two methods are shown to be able to compute a few extreme singular values of 2 50 × 2 50 matrix accurately in a few seconds on desktop computers. In this article, we propose a randomized tensor network algorithm for the computation of an SVD low-rank approximation of a matrix. As we will demonstrate through numerical experiments, our proposed algorithm manages to achieve the same accuracy up to 17 times faster than MALS-SVD. Moreover, our algorithm is able to retrieve accurate approximations for cases where both ALS-SVD and MALS-SVD fail to converge. More specifically, the main contributions of this article are twofold: 1. We present a fast algorithm that is able to convert a given sparse matrix into MPO form with machine precision accuracy. 2. We present a MPO version of the randomized SVD algorithm that can outperform the current state-of-the-art tensor algorithms [18] for computing lowrank matrix approximations of large-scale matrices. 1 This article is organized as follows. In Section 2, some basic tensor concepts and notation are explained. We introduce the notion of Matrix Product Operators in Section 3. Our newly proposed algorithm to convert a matrix into the MPO representation is presented in Section 4. In Section 5, we present our randomized algorithm to compute an SVD low-rank approximation of a given matrix in MPO form. Numerical experiments in Section 6 demonstrate both the fast matrix to MPO conversion as well as our randomized algorithm. We compare the performance of our matrix conversion algorithm with both the TT-SVD [22] and TT-cross [27] algorithms, while also comparing the performance of our randomized algorithm with the ALS-SVD and MALS-SVD algorithms [18]. Section 7 presents some conclusions together with an avenue of future research.
2. Tensor basics and notation. Tensors in this article are multi-dimensional arrays with entries either in the real or complex field. We denote scalars by italic letters, vectors by boldface italic letters, matrices by boldface capitalized italic letters and higher-order tensors by boldface calligraphic italic letters. The number of indices required to determine an entry of a tensor is called the order of the tensor. A dth order or d-way tensor is hence denoted A ∈ R I1×I2×···×I d . An index i k always satisfies 1 ≤ i k ≤ I k , where I k is called the dimension of that particular mode. We use the MATLAB array index notation to denote entries of tensors. Suppose that A is a 4way tensor with entries A(i 1 , i 2 , i 3 , i 4 ). Grouping indices together into multi-indices is one way of reshaping the tensor. For example, a 3-way tensor can now be formed from A by grouping the first two indices together. The entries of this 3-way tensor are then denoted by A([i 1 i 2 ], i 3 , i 4 ), where the multi-index [i 1 i 2 ] is easily converted into 1 MATLAB implementations of all algorithms are distributed under a GNU lesser general public license and can be freely downloaded from https://github.com/kbatseli/TNrSVD. a single index as i 1 + I 1 (i 2 − 1). Grouping the indices into [i 1 ] and [i 2 i 3 i 4 ] results in a Grouping indices together in order to change the order of a tensor is called reshaping and is an often used tensor operation. We adopt the MATLAB/Octave reshape operator "reshape(A, [I 1 , I 2 , · · · , I d ])", which reshapes the d-way tensor A into a tensor with dimensions I 1 × I 2 × · · · × I d . The total number of elements of A must be the same as I 1 × I 2 × · · · × I d . The mode-n matricization A (n) of a d-way tensor A maps the entry A(i 1 , i 2 , · · · , i d ) to the matrix element with row index i n and column index Example 1. We illustrate the reshaping operator on the 4 × 3 × 2 tensor A that contains all entries from 1 up to 24. Its mode-1 matricization is Another important reshaping of a tensor A is its vectorization, denoted vec(A) and obtained from grouping all indices into one multi-index.

Suppose we have two d-way tensors
Similarly, the outer product D = A • B of the d-way tensors A, B is a 2d-way tensor of dimensions I 1 × · · · × I d × J 1 × · · · × J d such that From equations (2.2) and (2.3) one can see that the Kronecker and outer products are interrelated through a reshaping and a permutation of the indices. A very convenient graphical representation of d-way tensors is shown in Figure 2.1. Tensors are here represented by circles and each 'leg' denotes a particular mode of the tensor. The order of the tensor is then easily determined by counting the number of legs. Since a scalar is a zeroth-order tensor, it is represented by a circle without any lines. One of the most important operations on tensors is the summation over indices, also called contraction of indices. For example, the following mode product [17] of a 3-way tensor A ∈ R I×I×I with a matrix U 1 ∈ R R×I and a vector u 3 ∈ R I 2 also illustrates a simple tensor network, which is a collection of tensors that are interconnected through contractions. The tensor network in Figure 2.2 has two legs, which indicates that the network represents a matrix. This article uses a very particular tensor network structure, the Matrix Product Operator structure.

Matrix Product Operators.
In this section, we give a brief introduction to the notion of MPOs. Simply put, an MPO is a linear chain of 4-way tensors that represents a matrix and was originally used to represent an operator acting on a multi-body quantum system. Since their introduction to the scientific community in 2010 [24], MPOs have found many other applications. We now discuss the MPO representation of a matrix through an illustrative example. Suppose that we have a matrix A of size I 1 I 2 I 3 I 4 × J 1 J 2 J 3 J 4 , as shown in Figure 3.1. This matrix A can be represented by an MPO of four 4-way tensors, where the first tensor has dimensions Similarly, the kth tensor in the MPO of Figure 3.1 hence has dimensions R k × I k × J k × R k+1 . We require that R 1 = R 5 = 1, which ensures that the contraction of this particular MPO results in a 8-way tensor with entries A(i 1 , j 1 , i 2 , j 2 , i 3 , j 3 , i 4 , j 4 ). This tensor can then be permuted and reshaped back into the original matrix with entries A([i 1 i 2 i 3 i 4 ], [j 1 j 2 j 3 j 4 ]). The dimensions R k , R k+1 of the connecting indices in an MPO are called the MPO-ranks and play a crucial role in the computational complexity of our developed algorithms. The MPO-ranks are called canonical if they attain their minimal value such that the MPO represents a given matrix exactly. A very special MPO is obtained when all MPO-ranks are unity. The contraction of a rank-one MPO corresponds with the outer product of the individual MPO-tensors. Indeed, suppose that the MPO-ranks in Figure 3.1 are all unity. The 4-way tensors of the MPO are then reduced to matrices such that we can write A (k) := A (k) (1, :, :, 1) (k = 1, . . . , 4) and which is exactly the outer product of the matrices A (1) , A (2) , A (3) with A (4) into a 8-way tensor. The relation between the Kronecker and outer products, together with the previous example leads to the following important theorem.
It is important to note that the order of the MPO-tensors is reversed with respect to the order of the factor matrices in the Kronecker product. This means that the last factor matrix A (1) in the Kronecker product of Theorem 3.1 is the first tensor in the corresponding MPO representation. Theorem 3.1 can also be written in terms of the matrix entries as As mentioned earlier, the MPO-ranks play a crucial role in the computational complexity of the algorithms. For this reason, only MPOs with small ranks are desired. An upper bound on the canonical MPO-rank R k for an MPO of d tensors for which R 1 = R d+1 = 1 is given by the following theorem.

in [27]) For any matrix
Proof. The upper bound on the canonical MPO-rank R k can be determined from contracting the first k − 1 tensors of the MPO together and reshape the result into the Similarly, we can contract the kth tensor of the MPO with all other remaining tensors and reshape the result into the R k × I k J k · · · I d J d matrix C. The canonical rank R k is now upper bounded by the matrix product BC due to rank(BC) ≤ min(rank(B), rank(C)).
These upper bounds are quite pessimistic and are attained for generic tensors. For example, a generic full-rank I 10 × I 10 matrix has an exact MPO representation with a canonical MPO-rank R 6 = I 10 . This implies that any MPO representation with R 6 < I 10 will consequently be an approximation of the original matrix. Theorem 3.2 therefore allows us to conclude that MPOs are only useful when either the canonical MPO-ranks are small or a low-rank MPO exists that approximates the underlying matrix sufficiently.
Two important MPO operations are addition and rounding, which are easily generalized from MPS addition and rounding. Indeed, by grouping the I k and J k indices together, one effectively transforms the MPO into an MPS such that MPS addition and rounding can be applied. The addition of two MPOs results in an MPO for which the corresponding ranks, except R 1 and R d+1 , are added. The rounding operation repeatedly uses the singular value decomposition (SVD) on each of the MPO-tensors going from left-to-right and right-to-left in order to truncate the ranks R k such that a specific relative error tolerance is satisfied. For more details on MPS addition and rounding we would like to refer the reader to [22, p. 2305] and [22, p. 2308], respectively. An alternative MPO rounding operation that does not require the SVD but relies on removing parallel vectors is described in [16]. This alternative rounding procedure can be computationally more efficient for MPOs that consist of sparse tensors.
4. Converting a sparse matrix into an MPO. The standard way to convert a matrix into MPO form is the TT-SVD algorithm [22, p. 2301], which relies on consecutive reshapings of the matrix from which an SVD needs to be computed. This procedure is not recommended for matrices in real applications for two reasons. First, application-specific matrices tend to be sparse and computing the SVD of a sparse matrix destroys the sparsity, which results in requiring more and often prohibitive storage. Second, real life matrices are typically so large that it is infeasible to compute their SVD. An alternative method to convert a matrix into an MPO is via cross approximation [27]. This method relies on heuristics to find subsets of indices for all modes of a tensor in order to approximate it. In practice, the cross approximation method can be very slow and not attain the desired accuracy. Our matrix to MPO conversion method relies on a partitioning of the sparse matrix such that it is easily written as the addition of rank-1 MPOs.

Algorithm derivation.
We derive our algorithm with the following illustrative example. Suppose we have a sparse matrix A ∈ R I×J with the following sparsity pattern Assume that each of the nonzero block matrices has dimensions I 1 × J 1 and that I = I 1 I 2 , J = J 1 J 2 such that the rows and columns of A are now indexed by [ respectively. The main idea of our method is to convert each nonzero block matrix into a rank-1 MPO and add them all together. Observe now that from which we conclude that the MPO representation of A is found from adding the unit-rank MPOs of each of the terms. Another important conclusion is that the MPO-rank for the particular MPO obtained from this algorithm is the total number of summations. The number of factors in the Kronecker product is not limited to two and depends on the matrix partitioning. Indeed, suppose we can partition A 14 further into then the first term of (4.1) becomes E 14 ⊗ E 21 ⊗ X 14 and likewise for the other terms.
A crucial element is that the matrix A is partitioned into block matrices of equal size, which is required for the addition of the MPOs. In general, for a given matrix A ∈ R I×J , we consider the partitioning of A determined by a Kronecker product of d matrices The algorithm to convert a sparse matrix into an MPO is given in pseudo-code in Algorithm 4.1.
Construct rank-1 MPO T with X and E ij matrices.
A ← A + T end for 4.2. Algorithm properties. Having derived our sparse matrix conversion algorithm, we now discuss some of its properties. Algorithm 4.1 has the following nice features, some of which we will address in more detail: • Except for A (1) , almost all of the MPO-tensors will be sparse.
• The user is completely free to decide on how to partition the matrix A, which determines the number of tensors in the resulting MPO. • The generalization of Algorithm 4.1 to construct an MPO representation of a given tensor is straightforward. • The maximal number of tensors in an MPO representation are easily deduced and given in Lemma 4.1. • The obtained MPO-rank for a particular partitioning of the matrix A is also easily deduced and given in Lemma 4.2. • A lower bound on the obtained MPO-rank for a fixed block size I 1 , J 1 is derived in Lemma 4.3. • As the dimensions of each of the MPO-tensors are known a priori, one can preallocate the required memory to store the tensors in advance. This allows a fast execution of Algorithm 4.1. The maximal number of tensors in an MPO representation of a matrix is determined by choosing a partitioning such that each block matrix of A becomes a single scalar entry and is given by the following lemma.
Lemma 4.1. Given a matrix A ∈ R I×J , suppose d I , d J are the number of factors in the prime factorizations of I, J, respectively. Then the maximal number of tensors in an MPO representation of A is max (d I , d J ) + 1.
The following simple example illustrates the maximal number of tensors max(d I , d J )+ 1 from Lemma 4.1.
The rank of the MPO obtained from Algorithm 4.1 for a particular partitioning of the matrix A is given by the following lemma.
Lemma 4.2. The MPO obtained from Algorithm 4.1 has a uniform MPO-rank equal to the total number of nonzero matrix blocks X as determined by the partitioning of A.
Lemma 4.2 follows trivially from the fact that ranks are added in MPO addition and all MPOs in Algorithm 4.1 are unit-rank. It is important to realize that the usage of Algorithm 4.1 is not limited to sparse matrices per se. One could apply Algorithm 4.1 to dense matrices but then the possible computational benefit of having to process only a few nonzero matrix blocks X is lost. It is also the case that the MPO-ranks can be reduced in almost all cases via a rounding procedure without the loss of any accuracy, as the upper bounds of Theorem 3.2 are usually exceeded. Partitioning the matrix A such that each term in Algorithm 4.1 corresponds with a single scalar entry sets the resulting MPO-rank to the total number of nonzero entries of A. This might be too high in practice. On the other hand, choosing any I k , J k too large results in a large MPO-tensor A (k) , which is also not desired. A strategy that can work particularly well is to use the Cuthill-Mckee algorithm [7] to permute A into a banded matrix with a small bandwidth. Grouping all nonzero entries together around the main diagonal also effectively reduces the number of nonzero block matrices and hence the total MPO-rank. A block size I 1 × J 1 can then be chosen such that the bandwidth is covered by a few blocks. Algorithm 4.1 can then be applied to the permuted matrix. Other permutations may reduce the maximal MPO-rank even further. We will discuss choosing the partitioning of A in more detail in Section 4.3. Algorithm 4.1 will construct an MPO with a uniform MPO-rank, which will exceed the upper bounds from Theorem 3.2 in almost all cases. One can use a rounding step to truncate the MPO-ranks without the loss of any accuracy after Algorithm 4.1 has finished. Alternatively, one can apply a rounding step on the intermediate result as soon as the MPO-rank reaches a certain threshold during the execution of the algorithm. The following example illustrates the necessity of the rounding step.
Example 4. Suppose we have three random 2 × 2 matrices A (1) , A (2) , A (3) such that By Theorem 3.1, the matrix A has a canonical unit-rank MPO representation where the first MPO-tensor is A (1) reshaped into a 1 × 2 × 2 × 1 tensor. Choosing I 1 = J 1 = I 2 = J 2 = I 3 = J 3 = 2 and applying Algorithm 4.1 results in an MPO with a uniform rank of 16. Applying a rounding step truncates each of these ranks down to unity.
For a fixed block size I 1 , J 1 one can determine a lower bound for the resulting MPO-rank in the following manner. Proof. Suppose that we found a permutation such that all z nonzero entries can be arranged into R block matrices of size I 1 × J 1 . It then trivially follows that R will be the MPO-rank since z = I 1 × J 1 × R.
In practice, it will be difficult, or in some cases impossible, to find a permutation such that all nonzero entries are nicely aligned into I 1 × J 1 block matrices. The R in Lemma 4.3 is therefore a lower bound.

Choosing a partition.
In this subsection we discuss choosing a partition of the matrix A prior to applying Algorithm 4.1. We suppose, without loss of generality, that the dimensions of A have prime factorizations I = I 1 I 2 · · · I d and J = J 1 J 2 · · · J d with an equal amount of d factors. The number of factors can always be made equal by appending ones. Ultimately, the goal is to obtain an MPO with "small" MPO-ranks. Although Algorithm 4.1 constructs an MPO with ranks that are likely to exceed the canonical values, these ranks can always be truncated through rounding. Theorem 3.2 can be used for choosing a partition that minimizes the upper bounds in the hope that the canonical values are even smaller. The key idea is that the upper bounds depend on the ordering of the prime factors. The following small example illustrates.
In light of the randomized SVD algorithm that is developed in Section 5 it will be necessary to order the prime factors in a descending sequence. In this way, the first MPO-tensor will have sufficiently large dimension in order to compute the desired low-rank approximation. Observe that if we use the descending ordering I 1 = 7, J 1 = 3, I 2 = 5, J 2 = 2, I 3 = 1, J 3 = 2 in Example 5, then the upper bounds are R 2 ≤ min(21, 20) = 20 and R 3 ≤ min(210, 2) = 2. Large matrices can have dimensions with a large number of prime factors. Choosing a partition with a large number of MPO-tensors usually results in a high number of nonzero block matrices X and therefore also in a large MPO-rank. The problem with such a large MPO-rank can be that it becomes infeasible to do the rounding step due to lack of sufficient memory. In this case one needs to reduce the number of MPO-tensors until the obtained MPOrank is small enough such that the rounding step can be performed. This way of choosing a partition will be demonstrated in more detail by means of a worked-out example in Section 6.1.

Tensor network randomized SVD.
5.1. The rSVD algorithm. Given a matrix A, which does not needs to be sparse, the rSVD computes a low-rank factorization U SV T where U , V are orthogonal matrices and S is a diagonal and nonnegative matrix. The prototypical rSVD algorithm [15, p. 227] is given as pseudocode in Algorithm 5.1. When a rank-K/2 approximation is desired, we compute a rank-K approximation after which only the first K/2 singular values and vectors are retained. This is called oversampling and for more details on this topic we refer the reader to [15, p. 240].It has been shown that a slow decay of the singular values of A results in a larger approximation error. The power iteration (AA T ) q tries to alleviate this problem by increasing the decay of the singular values while retaining the same left singular vectors of A. Common values for q are 1 or 2. Note that the computation of Y is sensitive to round-off errors and additional orthogonalization steps are required. For a large matrix A, it quickly becomes infeasible to compute orthogonal bases for Y or to compute the SVD of B. This is the main motivation for doing all steps of Algorithm 5.1 in an MPOform. We therefore assume that all matrices in Algorithm 5.1 can be represented by an MPO with relatively small MPO-ranks. For a matrix A ∈ R I1···I d ×J1···J d with d MPO-tensors A (i) ∈ R Ri×Ii×Ji×Ri+1 (i = 1, . . . , d), Algorithm 5.1 computes a rank-K factorization that consists of d MPO-tensors U (1) , . . . , U (d) and V (1) , . . . , V (d) and the K × K diagonal and nonnegative S matrix. We denote this MPO-version of the rSVD algorithm the tensor network randomized SVD (TNrSVD). The rSVD algorithm relies on multiplying the original matrix A with a random matrix O. Fortunately, it is possible to directly construct a random matrix into MPO form.
Lemma 5.1. A particular random J 1 J 2 · · · J d × K matrix O with rank(O) = K and K ≤ J 1 can be represented by a unit-rank MPO with the following random MPOtensors Proof. All MPO-ranks being equal to one implies that Theorem 3.1 applies. The random matrix O is constructed from the Kronecker product of d − 1 random column vectors O (i) ∈ R Ji , (i = 2, . . . , d) with the matrix O (1) ∈ R J1×K . The Kronecker product has the property that rank(O) = rank(O (d) ) · · · rank(O (1) ).
The fact that O (1) is a random matrix then ensures that rank(O) = K.
Probabilistic error bounds for Algorithm 5.1 are typically performed for random Gaussian matrices O [15, p. 273]. Another type of test matrices O are subsampled random Fourier transform matrices [15, p. 277]. The random matrix in MPO form from Theorem 5.1 will not be Gaussian, as the multiplication of Gaussian random variables is not Gaussian. This prevents the straightforward determination of error bounds for the MPO-implementation of Algorithm 5.1 that we propose. In spite of the lack of any probabilistic bounds on the error, all numerical experiments that we performed demonstrate that the orthogonal basis that we obtain for the range of A can capture the action of A sufficiently. Once the matrix A has been converted into an MPO using Algorithm 4.1 and a random MPO has been constructed using Theorem 5.1, what remains are matrix multiplications and computing low-rank QR and SVD factorizations. We will now explain how these steps can be done efficiently using MPOs.

Matrix multiplication. Matrix multiplication is quite straightforward. Suppose the matrices
flops. Corresponding MPO-ranks R i , S i and R i+1 , S i+1 are multiplied with one another, which necessitates a rounding step in order to reduce the dimensions of the resulting MPO-tensors. Note, however, that the random matrix O constructed via Lemma 5.1 has a unit-rank MPO, which implies that S 1 = S 2 = · · · = S d+1 = 1 such that the MPO corresponding with the matrix AO will retain the MPO-ranks of A.

Thin QR and economical SVD in MPO-form.
An orthogonal basis for the range of Y ∈ R I×K can be computed through a thin QR decomposition Y = Q R, where Q ∈ R I×K has orthogonal columns and R ∈ R K×K . The algorithm to compute a thin QR decomposition from a matrix in MPO-form is given in pseudo code in Algorithm 5.2. The thin QR is computed by an orthogonalization sweep from right-to-left, absorbing the R factor matrix into the preceding MPO-tensor. The main operations in the orthogonalization sweep are tensor reshaping and the matrix QR decomposition. The first MPO-tensor is orthogonalized in a slightly different way such that the K × K R matrix is obtained. The computational cost of Algorithm 5.2 is dominated by the QR computation of the first MPO-tensor, as it normally has the largest dimensions. Using Householder transformations to compute this QR decomposition costs approximately O(I 2 1 R 2 2 K) flops. The proof of the procedure can be The rSVD algorithm also requires an economical SVD computation of the K × J 1 · · · J d matrix B = W S V T , where both W , S are K × K matrices, W is orthogonal and S is diagonal and nonnegative. The matrix V is stored in MPO-form. Only a slight modification of Algorithm 5.2 is required to obtain the desired matrices. Indeed, the only difference with Algorithm 5.2 is that now the SVD of A 1 needs to be computed. From this SVD we obtain the desired W , S matrices and can reshape and permute the right singular vectors into the desired V (1) MPO-tensor. Again, the overall computational cost will be dominated by this SVD step, which costs approximately O(J 2 1 R 2 2 K) flops. The economical SVD of a matrix in MPO-form is given in pseudo-code in Algorithm 5.3. A graphical representation of the corresponding tensor network for a simple example of 4 MPO-tensors is depicted in Figure 5.3.
: Thin QR decomposition QR as a tensor network with Q ∈ R I1I2I3I4×K and R ∈ R K×K .
Both the thin QR and economical SVD of the matrix in MPO-form are computed for each tensor of the MPO separately, which reduces the computational complexity significantly. Unlike the ALS-SVD and MALS-SVD, no iterative sweeping over the different MPO-tensors is required.

Randomized subspace iteration.
The computation of the matrix Y = (AA T ) q AO is vulnerable to round-off errors and an additional orthogonalization step is required between each application of A and A T [15, p. 227]. Instead of computing Y and applying Algorithm 5.2, the randomized MPO-subspace iteration of Algorithm 5.4 is proposed. First, the random matrix O is multiplied onto A, after which a rounding step is performed to reduce the MPO-ranks. Algorithm 5.2 is then applied to obtain an orthogonal basis Q for the range of Y . One now proceeds with the multiplication A T Q, after which another rounding step and orthogonalization through Algorithm 5.2 are performed. These steps are repeated until the desired number of multiplications with A and A T have been done. The SVD-based rounding and orthogonalization steps can actually be integrated into one another. Indeed, one can apply a left-to-right rounding sweep first, followed by the right-to-left sweep of    and the bandwidth is reduced to 1962. We can factor 150102 as 2 × 3 × 3 × 31 × 269, which sets the maximal number of tensors in the MPO to 6. The total number of nonzero entries is 726674, which makes an MPO representation of 6 tensors infeasible as for this case all R k = 726674 and there is insufficient memory to store the MPO-tensors. Table 6.1 lists the number of MPO-cores d, the obtained MPOrank for both the original matrix and after applying the Cuthill-Mckee algorithm and the runtime for applying Algorithm 4.1 on the permuted matrix. Applying the permutation effectively reduces the MPO-rank approximately by half so we only consider the permuted matrix. First, we order the prime factors in a descending fashion 269, 31, 3, 3, 2, which would result in an MPO that consists of the following five ten- Due to the high MPO-rank, however, it is not possible to construct this MPO. We can now try to absorb the prime factor 2 into 269 and construct the corresponding MPO that consists of four tensors with A (1) ∈ R 1×538×538×1753 . This takes about 5 seconds. As Table 6.1 shows, increasing the dimensions of A (1) by absorbing it with more prime factors further reduces the MPO-rank and runtimes. Note that if MATLAB supported sparse tensors by default, then it would be possible to use Algorithm 4.1 for both the original and permuted matrix as all MPO-tensors are sparse. SVD-based rounding on the d = 3 MPO with a tolerance of 10 −10 reduces the MPO-rank from 347 down to 7 and takes about 61 seconds. Using the alternative parallel vector rounding from [16] truncates the MPO-rank also down to 7 in about 5 seconds.
6.2. Fast matrix-to-MPO conversion. In this experiment, Algorithm 4.1 is compared with the TT-SVD algorithm [24, p. 2135] and the TT-cross algorithm [27, p. 82], [28], both state-of-the-art methods for converting a matrix into an MPO. Both the TT-SVD and TT-cross implementations from the TT-Toolbox [23] were used. The TT-cross method was tested with the DMRG CROSS function and was run for 10 sweeps and with an accuracy of 10 −10 We use the 15838 × 15838 power simulation matrix from Liu Wenzhuo, EPRI, China, also in the UF Sparse Matrix Collection [8]. The prime factorization of 15838 consists of only two prime factors, 2 and 7919, and is not very suitable for a MPO conversion. We therefore append the matrix with zeros such that its dimensions are rounded to the nearest power of 2, as the SVD of the original matrix is easily recovered from the appended matrix. An MPO of 7 tensors is then constructed from the appended 16384 × 16384 matrix with dimensions I 1 = J 1 = 256, I 2 = · · · = I 7 = J 1 = · · · = J 7 = 2. Applying the Cuthill-Mckee algorithm reduces the maximal MPO-rank from 413 to 311 and is therefore not applied as the resulting rank decrease is not very significant. Table 6.2 lists the runtimes and relative errors when converting the appended matrix into the MPO format for the three considered methods. The relative errors are obtained by contracting the obtained MPO back into a matrixÂ and computing ||A −Â|| F /||A|| F . Both Algorithm 4.1 and the TT-SVD manage to obtain a result that is accurate up to machine precision. Although the maximal MPO-rank obtained with the TT-SVD algorithm is 22, which is an order of magnitude smaller than 413, our proposed algorithm is about 509 times faster than the TT-SVD algorithm. The TT-cross method fails to find a sufficiently accurate MPO and takes about the same amount of time as the TT-SVD algorithm.
Applying the TT-SVD and TT-cross method on the AMD-G2-circuit matrix was not possible due to insufficient memory.

Comparison with ALS-SVD and MALS-SVD.
The tensor networkmethod described in [18] uses the alternating least squares (ALS) and modified alternating least squares (MALS) methods to compute low-rank approximations of a given matrix in MPO-form. Three numerical experiments are considered in [18], of which two deal with finding a low-rank approximation to a given matrix. The first matrix that is considered is a rectangular submatrix of the Hilbert matrix. The Hilbert matrix H ∈ R 2 N ×2 N is a symmetric matrix with entries H(i, j) = (i + j − 1) −1 , i, j = 1, 2, . . . , 2 N . For this experiment, the submatrix A := H(:, 1 : 2 N −1 ) is considered with 10 ≤ N ≤ 50. Following [18], the corresponding MPO is constructed using the FUNCRS2 function of the TT-Toolbox [23], which applies a TT-cross approximation method using a functional description of the matrix entries, and consists of N MPOtensors. The obtained MPO approximates the Hilbert matrix with a relative error of 10 −11 . Maximal MPO-ranks for all values of N were bounded between 18 and 24. A tolerance for the relative residual of 10 −8 was set for the computation of rank-16 approximations with both the ALS and MALS algorithms. In order to be able to apply Algorithm 5.1, we first need to make sure that the MPO-tensor Y (1) of the matrix Y = AO has dimensions 1 × I 1 × K × R 2 with I 1 ≥ K. Since a rank-16 approximation is desired, this means that K = 32. By contracting the first 5 MPOtensors A (1) , A (2) , A (3) , A (4) , A (5) into a tensor with dimensions 1 × 32 × 32 × R 6 , we obtain a new MPO of N − 5 + 1 tensors that satisfies the I 1 ≥ K condition. A tolerance of 10 −9 was used in the rounding procedure and q was set to 2, 3, 4 for N = 10, . . . , 30, N = 35, 40 and N = 45, 50, respectively. Figure 6.3 shows the runtimes of the ALS, MALS and TNrSVD method as a function of N . All computed rank-16 approximations obtained from the ALS and MALS methods satisfied ||A − U SV T || F /||A|| F ≤ 10 −8 , as reported by the respective methods. As mentioned in [18], computing the norm of the residual in MPO-form can be computationally challenging when the MPO-ranks are high. For this reason, we compared the obtained singular values from the ALS method with the singular values from the TNrSVD method to ensure the relative residuals were below 10 −8 . The MALS method solves larger optimization problems than the ALS method at each iteration and is therefore considerably slower. The TNrSVD method is up to 6 times faster than the ALS method and 13 times faster than the MALS method for this particular example. Using a standard matrix implementation of Algorithm 5.1 we could compute low-rank approximations only for the N = 10 and N = 15 cases, with respective runtimes of 0.04 and 9.79 seconds. From this we can conclude that all three tensor-based methods outperform the standard matrix implementation of Algorithm 5.1 when N ≥ 15 for this particular example. The second matrix that is considered is one with 50 prescribed singular values 0.5 k , (k = 0, . . . , 49) and random left and right singular vectors U ∈ R 2 N ×50 , V ∈ R 2 N ×50 . As with the Hilbert matrices, N ranges from 10 up to 50 and equals the number of tensors in the MPO. The maximal MPO-rank of all constructed MPOs was 25. The orthogonal U , V matrices were generated in MPO-form using the TT RAND function of the TT-Toolbox. The MPO-representation of the matrix was then obtained by computing U SV T in MPO-form. A tolerance of 10 −6 was set for the computation of rank-50 approximations with both the ALS and MALS algorithms such that all approximations satisfy ||S −Ŝ|| F /||S|| F ≤ 10 −6 , whereŜ denotes the diagonal matrix obtained from either the ALS or MALS method. One sweep of the ALS and MALS algorithms sufficed to obtain the desired accuracy, except for the cases N = 45 and N = 50. For these two cases neither ALS nor MALS was able to converge to the desired accuracy. The exponent q was set to 1 for all runs of the TNrSVD algorithm and the rounding tolerances were set to 10 −5 , 10 −6 , 10 −8 , 10 −8 , 10 −9 , 10 −10 , 10 −11 , 10 −12 , 10 −13 for N = 10, 15, 20, . . . , 50, respectively. This ensured that the result of the TNrSVD method had a relative error ||S −Ŝ|| F /||S|| F on the estimated 50 dominant singular values below 10 −6 . Computing a rank-50 approximation implies that K = 100 and the first 7 MPO tensors need to be contracted prior to running the TNrSVD algorithm. These contractions result in a new MPO where the first tensor has dimensions 1 × 128 × 128 × R 8 , such that I 1 = 128 ≥ K = 100 is satisfied. Figure 6.4 shows the runtimes of the ALS, MALS and TNrSVD method as a function of N . Just like with the Hilbert matrices, the MALS algorithm takes considerately longer to finish one sweep. The TNrSVD algorith is up to 5 times faster than ALS and 17 times faster than MALS for this particular example. Using a standard matrix implementation of Algorithm 5.1 we could compute low-rank approximations only for the N = 10 and N = 15 cases, with respective runtimes of 0.04 and 7.59 seconds. For N = 15, the runtimes for the ALS, MALS and TNrSVD methods were 7.63, 19.5 and 2.13 seconds, respectively. 7. Conclusion. We have proposed a new algorithm to convert a sparse matrix into an MPO form and a new randomized algorithm to compute a low-rank approximation of a matrix in MPO form. Our matrix to MPO conversion algorithm is able to generate MPO representations of a given matrix with machine precision accuracy up to 509 times faster than the standard TT-SVD algorithm. Compared with the state-of-the-art ALS-SVD and MALS-SVD algorithms, our TNrSVD is able to find low-rank approximations with the same accuracy up to 6 and 17 times faster, respectively. Future work includes the investigation of finding permutations, other than the Cuthill-Mckee permutation, that can reduce the MPO-rank of a given matrix.