We study the limitations of matrix-free algorithms that access a matrix \(A\) only through forward matrix-vector products (matvecs) \(x \mapsto Ax\), without access to the transpose \(A^\top\) or its action. This setting arises naturally in operator learning, inverse problems, and matrix-free PDE solvers, where adjoint evaluations may be unavailable or prohibitively expensive. We show that the lack of transpose access creates severe and sometimes insurmountable theoretical barriers. For Krylov methods, we prove that the sequence of projected operator norms produced by Arnoldi iteration can follow any prescribed nondecreasing curve, showing that forward matvecs alone provide essentially no reliable information about the spectral norm. For several core problems, including least squares, norm estimation, column subset selection, and local maximum volume, we establish non-identifiability results; distinct matrices can generate identical forward-query transcripts while having fundamentally different solutions. We also prove quantitative lower bounds on the number of forward matvecs required for approximation tasks. In particular, any algorithm that computes a near-optimal rank-\(k\) approximation must use at least \(n\) queries, and estimating the Frobenius norm to relative accuracy \(\varepsilon\) requires \(\Omega(\varepsilon^{-2})\) queries when \(n\) is sufficiently large, matching the complexity of Hutchinson-type estimators up to constants. Although some problems remain solvable without transpose access, the transpose-free setting is fundamentally more limited in both identifiability and efficiency.

1 Introduction↩︎

In modern numerical linear algebra, it is common to treat a matrix \(A\) as a black box that can be queried through matrix-vector products. This matrix-free viewpoint underpins Krylov subspace methods, randomized linear algebra, and large-scale scientific computing more broadly, especially when forming or storing \(A\) explicitly is infeasible. In many applications, however, access to the transpose or adjoint operator is unavailable, inaccurate, or significantly more expensive than forward evaluations. This occurs in matrix-free discretizations of differential equations [1], operator learning from input-output data [2], [3], matrix recovery [4], [5], and inverse problems in which the adjoint corresponds to a different physical model. These examples raise a basic question: what can be computed about a matrix when one has access only to the map \(x \mapsto Ax\), but not to \(y \mapsto A^\top y\)?

We study this question in a query model. The matrix \(A \in \mathbb{R}^{m \times n}\) is accessed through an oracle that returns \(Ax\) for any query vector \(x \in \mathbb{R}^n\). An algorithm may choose queries \(x_1,\dots,x_q\) adaptively, with each \(x_i\) depending on previous queries and responses and possibly on internal randomness. After observing the transcript \((x_1,Ax_1),\dots,(x_q,Ax_q)\), the algorithm outputs an estimate of a target quantity, such as a norm, a singular vector, or a low-rank approximation. We call an algorithm transpose-free (TF) if it accesses \(A\) only through such forward queries and has no access to \(A^\top\) or to the entries of \(A\) itself. Our main object of study is the number of forward queries needed to solve a given problem to prescribed accuracy and success probability. This setup is intentionally minimal: we do not restrict the algorithmic paradigm, the amount of internal computation, or the adaptivity of the queries. The only resource we measure is access to the forward oracle.

This model lets us separate two distinct issues. The first is identifiability: whether an efficient set of forward queries uniquely determines the quantity of interest at all. The second is query complexity: how many forward matvecs are required to solve a problem to a given accuracy. This distinction turns out to be useful throughout the paper. Many problems fail at the level of identifiability; the same forward information is consistent with multiple incompatible answers arising from different, indistinguishable matrices. For some problems, we prove a query complexity lower bound showing that regardless of the transcript, a problem may be solvable with only forward matvecs, but only at a query cost that is substantially larger than what is possible when both \(A\) and \(A^\top\) are available.⁵ Our main message is that TF access is fundamentally limited in what can be computed efficiently compared to the standard matrix-free access setting with both \(A\) and \(A^\top\) available.

Our results fall into three categories. First, we prove non-identifiability results for several problems, including least squares, maximum volume, and aspects of spectral estimation. Given a ground-truth matrix \(A\) and a transcript, we construct an adversarial matrix \(B\) with an identical transcript, but yielding a very different solution. These are strong impossibility results in the sense that no amount of computation can recover information that is absent from the transcript. Second, we prove lower bounds on TF query complexity, independent of the transcript. In particular, a near-optimal rank-\(k\) approximation of an \(n \times n\) matrix requires \(n\) forward queries, and estimating the Frobenius norm to relative accuracy \(\varepsilon\) requires \(\Omega(\varepsilon^{-2})\) forward queries. Third, we show that TF Krylov information can be highly misleading. The projected norms produced by Arnoldi can follow any prescribed nondecreasing curve, so these quantities do not reliably track the spectral norm. These three types of results do not immediately imply one another, and indeed sometimes yield very different lower bounds. Thus, they each provide different qualitative descriptions of the limitations of the TF setting.

TF computation has a long history in numerical linear algebra, especially in the development of Krylov methods. Much of the existing literature focuses on algorithm design under one-sided access. Our goal here is complementary. Instead of proposing yet another TF algorithm, we ask what can and cannot be inferred from one-sided information in the first place, and how expensive successful inference must be when it is possible.

The rest of the paper is organized as follows. In 2, we review some existing algorithms and theoretical results for TF linear algebra. In 3, we describe the techniques used to prove our original results in the TF setting. In 4, we study TF approximation of the spectral norm, showing that using Krylov queries, projected norms can follow arbitrary nondecreasing curves in 4.1 and that the solution is not identifiable unless one uses \(n\) queries in 4.2. In 5, we prove a sharp lower bound on the number of queries needed for near-optimal low-rank approximation and show that the set of columns spanning the range of the matrix is not identifiable unless one uses \(n\) queries. In 6, we establish non-identifiability for local maximum volume. In 7, we study Schatten-norm estimation, and prove a query complexity lower bound for estimating the Frobenius norm, resulting in a rank-one corollary for all Schatten norms. In 8, we prove a non-identifiability result for overdetermined least squares problems with orthonormal design matrices. We conclude with open problems in 9.

Notation: Unless otherwise specified, \(\|\cdot\|\) denotes the spectral norm for matrices and the Euclidean norm for vectors. The Frobenius norm is denoted by \(\|\cdot\|_F\). The matrix \(A\) is \(m \times n\) or \(n \times n\), depending on the problem. We denote by \(A_k\) the best rank-\(k\) approximation to \(A\) given by the Eckart–Young theorem. The matrix trace is denoted by \(\operatorname{tr}(A)\). The identity matrix of appropriate dimension is denoted by \(I\), and its \(i\)th column by \(e_i\).

1.1 Main contributions↩︎

We now describe our main technical contributions. We summarize our theoretical contributions, along with existing results, in 1.

1. TF Krylov information can be arbitrarily misleading. In 2, we show that for any prescribed positive nondecreasing sequence \(0 < \sigma_1 \le \cdots \le \sigma_n\), there exists a matrix \(A\) for which the Arnoldi projected norms satisfy \(\|H_{k+1,k}\| = \sigma_k\) for \(1 \le k \le n-1\), while \(\|A\| = \sigma_n\). Thus, the projected norms seen by a TF Krylov method may not reveal any meaningful information about the spectral norm. In 3, we strengthen this by showing that the same pathology can occur even after fixing the eigenvalues of \(A\).

2. Near-optimal low-rank approximation is matvec-expensive without the transpose. 5 shows that any randomized, possibly adaptive TF algorithm that returns a \((1+\varepsilon)\)-optimal rank-\(k\) approximation for every \(n \times n\) input matrix must use \(n\) forward matvecs. In particular, one cannot obtain a genuinely matrix-free low-rank approximation algorithm from forward products alone.

3. Norm estimation exhibits sharp TF lower bounds. For Frobenius norm estimation, 9 proves that relative-error estimation within \(\varepsilon\) requires \(\Omega(\varepsilon^{-2})\) forward queries, for sufficiently large \(n\), matching Hutchinson-type methods up to constants and ruling out Hutch++-type \(\mathcal{O}(\varepsilon^{-1})\) improvements in the TF model. In particular, Hutchinson’s estimator for Frobenius norm estimation is optimal among TF methods. As an immediate consequence, 1 gives the same lower bound for estimating any Schatten norm, as all Schatten norms coincide on rank-one matrices.

4. Several basic problems are not identifiable from one-sided transcripts. In 6, we show that a set of columns spanning the range of the matrix is not identifiable without \(n\) forward matvecs. In 8, we prove that a solution to the local maximum volume problem for an \(m \times n\) matrix is also not determined by fewer than \(n\) forward matvecs. In 11, we prove a non-identifiability result for overdetermined least squares problems with orthonormal matrices, even when the adversarial matrix has condition number 1. Finally, a simple adversarial construction in 4 establishes non-identifiability for norm estimation.

5. A review of the TF literature. While TF methods have existed for a long time, the present work provides a unified perspective on this literature. Beginning with the historical origins and motivations for the TF setting, we discuss methods for large-scale linear systems, ill-posed problems, and model reduction in 2. We conclude by describing operator learning, as well as the closely related “unmatched transpose” setting.

2 A brief review of transpose-free methods↩︎

Many classical algorithms for nonsymmetric problems, such as bi-Lanczos-based methods and their descendants, rely on access to both \(A\) and \(A^\top\). By contrast, Arnoldi-based methods such as GMRES [12] use only forward products. A substantial literature from the late 1980s and 1990s sought to remove transpose access while retaining short recurrences, leading to methods such as CGS [13], BiCGSTAB [14], and TFQMR [15], as well as TF reformulations based on polynomial and squaring techniques [1]. The broader literature also contains “anything is possible” pathology results for Krylov methods, most famously for GMRES [11]; see also [16]–[19]. For spectral norm estimation, however, standard algorithms typically exploit the transpose, e.g., via Lanczos bidiagonalization or Krylov–Schur methods applied to \(A^\top \! A\) [20], [21], and even rough estimators often use iterations with both \(A\) and \(A^\top\). More recently, the TF setting has appeared in operator learning and data-defined operators [2], [3]. In this section, we give a short review of TF algorithms and theory in numerical linear algebra.

2.1 Motivations for TF methods↩︎

Traditionally, there are three contexts where TF methods are especially useful [1]. First, for the solution of nonlinear equations \(F(x) = 0\), with \(F: \mathbb{R}^n \to \mathbb{R}^n\), the action of the Jacobian \(x \mapsto DF(a) \, x\) is often used in inner iterations for the solution of a large-scale linear system, as in Newton–Krylov methods. This action is then approximated by a difference quotient \(\delta^{-1} \, (F(a+\delta \, x)-F(a))\) for a small \(\delta\) (typically, \(10^{-8}\)), while there is no similar expression for the transpose. Second, the action of the forward operator may be implemented in a routine or software package, such as in integral equations [1] or computed tomography (CT) [22]. Finally, TF methods are natural for popular matrix storage types (e.g., compressed row format) which enable efficient matrix-vector products. The transpose may be much more expensive, or require twice the storage (a compressed column format).

2.2 Linear systems↩︎

Consider the solution of a nonsingular linear system \(Ax = b\), where \(A \in \mathbb{R}^{n \times n}\). A variety of Krylov methods (see, e.g., [23] for an overview) search for an approximate solution in the Krylov spaces \(\mathcal{K}_k(A,b) := {\rm span}\{b, Ab, \dots, A^{k-1}b\}\). Some methods are explicitly designed for symmetric \(A\), such as CG, MINRES, and SYMMLQ. Some methods for nonsymmetric \(A\), such as GMRES and FOM, do not need the action of the transpose, at the price of storing basis vectors and (re)orthogonalization against those vectors. Two-sided methods, like BiCG, use the transpose. Short recurrences enable reduction in both the storage of basis vectors and (re)orthogonalization costs. To exploit a shadow Krylov space \(\mathcal{K}_k(A^\top, c)\) for a shadow vector \(c\), these methods use vectors of the form \(\psi(A^\top) \, c\), where \(\psi\) is a degree-\(k\) polynomial. However, the convergence of BiCG is typically quite irregular and non-monotonic.

CGS, BiCGSTAB, and the TFQMR methods seek to improve on the BiCG approach, first by stabilizing the convergence. These methods are also suited to avoiding transposes. The essential trick is to access these vectors only through inner products, as \((\psi(A^\top) \, c, \;\varphi(A) \, b) = (c, \;\psi(A) \, \varphi(A) \, b)\). Rather than working with a residual and its conjugate (as in BiCG), CGS uses a squared residual, which also avoids the transpose. TFQMR is related to CGS and has the advantage of avoiding irregular convergence behavior by quasi-minimizing a residual norm. Both BiCGSTAB and TFQMR are widely used methods in science and engineering.⁶

Various other methods following a similar strategy have been proposed. BiCGSTAB(\(\ell\)) [24] uses polynomials of degree \(\ell\) instead of one for the shadow residual to get an even smoother convergence. A lesser-known variant is QMRCGSTAB [25], which is a QMR type adaptation of BiCGSTAB, in the same way that TFQMR is of CGS. Some other variations are discussed in [26]. A relatively recent TF linear solver, based on an idea of the 1980s, is IDR(\(s\)) [27]. On two-sided methods, it has been shown in [28] that in a certain class of methods exploiting a three-term recurrence, none minimizes the residual norm; see also [29]. A review on developments of Krylov methods in the 21st century can be found in [30].

TF methods have been proposed for several generalized problems. For the block case \(AX = B\), where \(B\) comprises many right-hand sides, [31] presents a hybrid TF block method, which is claimed to use less memory and perform better than a TF block GMRES method. The TFQMR method can also be elegantly used for families of shifted linear systems [32].

Finally, the computation of matrix functions, \(x = f(A) \, b\), may be viewed as a generalization of linear systems, which are obtained by \(f(A) = A^{-1}\). Several aspects of TF subspace extraction methods are compared with those of two-sided approaches (using \(A^\top\)) in [33].

In the realm of TF lower bounds for linear systems, [10] recently established a query complexity lower bound of \(n\) forward queries. In the present work, we establish an identifiability lower bound of \(n\) forward queries for orthonormal least squares problems.

2.3 Linear ill-posed problems↩︎

Linear ill-posed problems are of the form \(Ax \approx b\), with \(A \in \mathbb{R}^{m \times n}\), where all three situations \(m=n\), \(m > n\), and \(m<n\) occur frequently in practice, and \(A\) is (near) singular. To get sensible approximate solutions, it is common to add a regularization term; a popular alternative is to use an iterative method as regularization, by stopping early enough in the process. The iteration number \(k\) then acts as a regularization parameter.

For square problems, both the TF method GMRES and the non-TF approach LSQR are popular. There is no single best method for all problems. However, as reported by experts and practitioners, the method of choice for a wide range of problems is LSQR [34]. One heuristic reason for this is that small singular values play an important role for these problems and are approximated well by LSQR Krylov spaces of the form \(\mathcal{K}_k(A^\top \! A, \, A^\top b)\), while GMRES Krylov spaces \(\mathcal{K}_k(A, b)\) by nature target eigenvalues. Indeed, the fact that the default method, LSQR, queries \(A^\top\) is in line with the general message of this paper.

For the overdetermined case (\(m > n\)), only methods involving Krylov spaces with \(A^\top\) (such as LSQR, LSMR, CGLS, CGNR, CGNE, and SYMMLQ), seem to be natural options. An alternative to iterative methods as a regularizer is to add a regularization term involving a 2-norm, 1-norm, or other \(p\)-norm. For these approaches, flexible Krylov methods are popular, including a variable right preconditioner to obtain the desired norm regularization. The main underlying mechanisms are flexible Arnoldi (which is TF) and flexible Golub–Kahan (Lanczos) bidiagonalization (exploiting \(A^\top\)); see, e.g., [35].

2.4 Model reduction of LTI dynamical systems↩︎

A standard form of a (single-input, single-output) linear time-invariant (LTI) dynamical system is \(\dot{x}(t) = A \, x(t) + b \, u(t)\), \(y(t) = c^\top x(t)\). Here \(b\), \(c \in \mathbb{R}^n\), \(u: \mathbb{R}\to \mathbb{R}\) is the input function, and \(y: \mathbb{R}\to \mathbb{R}\) is the output function. For sufficiently large \(|s|\) so that the series converges, the associated transfer function has the expansion \[\label{eq:mor1} H(s) = c^\top (s\,I-A)^{-1} \, b = s^{-1} \, (c^\top \,b) + s^{-2} \, (c^\top \!A\,b) + s^{-3} \, (c^\top \!A^2\,b) + \cdots\tag{1}\] Here, \(s \in i \, \mathbb{R}\) indicates the frequency. The terms \(s^{-j} \, c^\top \!A^{j-1}b\) are the moments. Since \(n\) is often large, dimension reduction is desirable; see, e.g., [36] for an overview of the contents of this subsection. A breakthrough of the 1990s [37] shows that an accurate reduction of the transfer function can be obtained by projecting onto Krylov spaces, rather than explicitly computing moments. One default technique is to use two-sided projections, including the transpose. Let the columns of \(V\) and \(W\) span \(\mathcal{K}_k(A, b)\) and \(\mathcal{K}_k(A^\top, c)\), respectively. For simplicity, we assume that the bases are biorthogonal (\(W^\top V = I\)); minor adaptations are needed in case both \(V\) and \(W\) have orthonormal columns. Then, the reduced transfer function \[\widehat H(s) = (c^\top V) \, (s\,I-W^\top\!AV)^{-1} \, W^\top b\] matches the first \(2k\) moments in the expansion of the transfer function.

An alternative to this two-sided approach is a TF one-sided projection, using only \(V\) with orthonormal columns (so \(V^\top V = I\)). The corresponding approximation \[\label{eq:mor2} \widetilde{H}(s) = (c^\top V) \, (s\,I-V^\top\!AV)^{-1} \, V^\top b\tag{2}\] matches only \(k\) moments. This usually means that the approximation is slightly worse, especially for the higher frequencies. However, such a TF approach may have other advantages over a two-sided scheme. Most notably, a one-sided projection 2 preserves stability if the field of values \(\mathcal{F}(A) := \{ \, x^\top Ax \, : \, x \in \mathbb{R}^n, \, \|x\| = 1 \, \}\) is located in the open left-half plane, since this implies that this also holds for the spectrum of \(V^\top \! AV\) for any \(V\) with orthonormal columns. In this case, the eigenvalues of \(V^\top\!AV\) are necessarily in the left-half plane.

Although 1 represents a Taylor expansion at \(s=\infty\), matching the first couple of moments by a projection often also yields good fits for smaller frequency values. However, for better fits near a given frequency \(s_0\), a common approach is to consider the shifted transfer function \(c^T ((s-s_0)\,I - (A-s_0 \, I))^{-1} \, b\) instead, leading to moments of the form \(c^T (A-s_0 I)^{-j} \, b\).

For generalized linear systems \(E \, \dot{x}(t) = A\,x(t) + b \, u(t)\), with \(E\) nonsingular, the associated transfer function is of the form \[H(s) = c^\top (s\,E-A)^{-1} \, b = c^\top (s\,I-E^{-1} A)^{-1} \, E^{-1}b,\] and there are similar one-sided and two-sided methods available as for the standard case.

2.5 Operator learning↩︎

Increasingly, scientists and engineers use deep learning methods to build reduced order models of dynamical systems and solution operators of PDEs. Neural operators [3], [38] are infinite-dimensional analogues of neural networks used to approximate these operators using data generated from physical experiments or numerical simulations. For linear solution operators, a natural idealized model accesses both forward queries \(f \mapsto \mathcal{A}f\) and adjoint queries \(g\mapsto \mathcal{A}^\ast g\). In most practical settings, however, only \(f \mapsto \mathcal{A}f\) is available.

If \(\mathcal{A}\) is the solution operator of a linear PDE, then access to the adjoint \(\mathcal{A}^\ast\) is usually unrealistic as the adjoint is the solution operator to a different PDE. Thus, the goal of recent work [2] is to build a theoretical framework for adjoint-free operator learning. The central contribution of this work is an adjoint-free operator learning algorithm whose convergence relies on prior knowledge about \({\rm range}(\mathcal{A}^\ast)\) encoded in a prior, self-adjoint operator \(L\). In the case of elliptic PDEs, elliptic regularity is used to obtain an algebraic adjoint-free convergence rate [2]. For more general classes of operators, larger sample complexity guarantees that depend on the smoothness of the operator have been derived [39]–[41].

Thus, while adjoint-free operator learning may be impossible or highly inefficient for general operators, operators arising from PDEs or dynamical systems may have structure that enables adjoint-free operator learning. With adjoint access, existing convergence guarantees for learning elliptic PDEs are exponential, rather than algebraic [42]. It is an open question whether one can close this theoretical gap between the adjoint and adjoint-free settings.

Motivated by the operator learning setting, the analogous problem of matrix recovery from matrix-vector products has also been studied in many works [4], [43]–[45]. In [43], a theoretical gap is established between the query complexity of matrix recovery using one-sided versus two-sided queries. To find a near-optimal approximation within a finite family \(\mathcal{F}\) to a general matrix, they establish a quadratic improvement in query complexity over the TF setting. The algorithm uses sketches on one side to identify the fraction of “bad” approximations that sketches on the other side can eliminate from \(\mathcal{F}\). This barrier may extend to other matrix recovery problems.

2.6 Unmatched transpose↩︎

The previous subsections focus on TF methods for various problems. There has also been considerable attention to problems involving an operator with an unmatched transpose, or adjoint mismatch. This usually arises from the fact that the forward and transpose operators correspond to separate and extensive software implementations, where the resulting actions are not close to transposes of one another. A typical motivation comes from CT, see, e.g., [46], where \(A\) corresponds to the projection and \(A^\top\) to the backprojection. An unmatched transpose may be viewed as having inexact information about \(A\)’s row space.

Consider \(A\) with unmatched transpose \(B\). Typically, \(B\) is a “quite inaccurate transpose of \(A\),” i.e., \(BA\) is a nonsymmetric operator with, besides positive eigenvalues, real but negative eigenvalues, complex eigenvalues with positive real part, and also complex eigenvalues with negative real part; see, e.g., [47]. The operator \(BA\) gives rise to interesting recent developments. In [22], [48], column- and row-oriented iterative methods are studied. BA-GMRES is the topic of [49]. For a least squares problem \(Ax \approx b\), it is shown in [50] that “not-normal equations” of the form \(BAx = Bb\) may have conditioning advantages compared to the usual normal equations. In [6], the authors address the question of how to quantify the norm of the mismatch. A proximal gradient method for an ill-posed problem with a nonsmooth regularization term involving an unmatched transpose is discussed in [51].

3 Query model and identifiability↩︎

We now introduce the TF query model and the basic mechanisms that underlie the impossibility and lower-bound results in the rest of the paper. We use an oracle model closely related to standard models in matrix sketching and matrix recovery [4], [5], [52]. An unknown matrix \(A \in \mathbb{R}^{m\times n}\) is accessed only through the operation \[x \longmapsto Ax, \qquad x \in \mathbb{R}^n.\] At round \(i\), the algorithm chooses a query vector \(x_i\), possibly adaptively and using internal randomness, and receives the response \(Ax_i\). After \(q\) rounds, it has observed the transcript \[(x_1,Ax_1), \dotsc, (x_q,Ax_q),\] and must output the object of interest. The central structural fact is that the transcript depends only on the action of \(A\) on the query space \[\mathcal{X}_q = \mathrm{span}\{x_1,\dots,x_q\}.\] The action of \(A\) on anything outside of \(\mathcal{X}_q\) is invisible to the algorithm. This leads to two complementary notions that recur throughout the paper. First, a solution to a problem may be non-identifiable from the transcript in the sense that for every matrix \(A\), there exists a matrix \(B\) such that \[Ax = Bx \qquad \text{for all } x \in \mathcal{X}_q,\] even though the solutions to the problems involving \(A\) and \(B\) are significantly different. In this type of result, the construction of the adversary \(B\) depends on the transcript. As a consequence, no deterministic TF algorithm can solve the problem from the transcript alone. Second, while a randomized algorithm may avoid the adversary with high probability, it may still suffer from high query complexity; one needs many queries before the transcript contains enough information to guarantee success.

Importantly, a lower bound on query complexity does not imply a lower bound on the number of queries needed for identifiability, nor vice versa. For example, in 4, we show that for any matrix \(A\), one needs \(n\) queries with \(A\) for the Frobenius norm of \(A\) to be uniquely determined. However, Hutchinson’s estimator applied to \(A^\top A\) needs only \(\mathcal{O}(\varepsilon^{-2})\) queries to achieve an \(\varepsilon\)-approximation. This is because given a query space \(\mathcal{X}_q\), a non-identifiability argument constructs an adversary using \(\mathcal{X}_q^\perp\). Randomized algorithms can circumvent this obstruction because \(\mathcal{X}_q\) itself is random, and a fixed adversary is unlikely to remain indistinguishable from \(A\) across all random query realizations. Conversely, query complexity lower bounds obtained via Yao’s minimax principle do not explicitly characterize the adversarial matrices associated with \(\mathcal{X}_q\), whereas non-identifiability arguments directly exhibit such matrices, providing geometric insight into the problem. Thus, these types of results are qualitatively different measures of hardness and formally incomparable.

3.1 Yao’s minimax principle and hard distributions↩︎

Our quantitative lower bounds on query complexity are proved using Yao’s minimax principle [53]. In the present setting, the principle says that to lower-bound the query complexity of randomized, possibly adaptive TF algorithms, it suffices to find a hard distribution of matrices and analyze the average performance of all deterministic, adaptive algorithms over that distribution. Formally, if \(\mathcal{A}\) denotes a class of deterministic algorithms, \(\mathcal{R}\) the associated randomized algorithms, and \(c(\mathrm{alg},A)\) a cost function, then \[\inf_{R \, \in \, \mathcal{R}} \sup_{A} \;\mathbb{E}[c(R, A)] \;\ge\; \sup_{\mathcal{D}} \inf_{\mathrm{alg} \, \in \, \mathcal{A}} \;\mathbb{E}_{A \sim \mathcal{D}}[c(\mathrm{alg}, A)].\] Thus, once a hard distribution \(\mathcal{D}\) is fixed, one may treat the algorithm as deterministic and focus entirely on what the transcript can and cannot reveal.

In Sections 5 and 7, we use two kinds of cost functions. For low-rank approximation, we take \(c(\mathrm{alg},A)\) to be the failure indicator, so that its expectation is the failure probability. For Schatten-norm estimation, we instead take \(c(\mathrm{alg},A)\) to be the estimation error. The lower-bound arguments then follow a common pattern: construct a hard distribution, consider any deterministic algorithm with \(q\) queries, and show that after observing the full transcript, there remains either ambiguity in the correct answer or enough residual uncertainty to force a nontrivial error. In this paper, we appeal to Yao’s minimax principle in proving lower bounds on the TF query complexity of low-rank approximation and Frobenius/Schatten norm estimation.

3.2 The mechanism behind identifiability↩︎

Several of our negative results are driven by the same geometric mechanism. We say that the solution to the problem is non-identifiable if for all \(A\), there exists a distinct matrix \(B\) which satisfies the same matvecs as \(A\), i.e., \(AX = BX\), but the solutions to the problems with \(A\) and \(B\) are very different. The geometric insight is that one can “hide” directions orthogonal to \(X\) in the matrix \(B\) to produce a very different answer. This is formalized in the following proposition.

Thus, determining if the solution of a linear algebra problem is identifiable from the transcript is equivalent to finding a \(Z= B\,Q_2\) that results in a drastically different solution. For example, one may construct \(Z\) so as to alter \(A\)’s columns lying outside of \(\mathcal{X}_q\), change the least squares solution by acting in a direction orthogonal to \(\mathcal{X}_q\), or plant additional norms or ranks in a direction orthogonal to \(\mathcal{X}_q\). This is the basic source of the non-identifiability results in later sections; the algorithm learns only a partial action of the matrix, and different global matrices can agree on that partial action.

4 Spectral norm estimation↩︎

In this section, we consider the problem of matrix-free spectral norm estimation. With access to the transpose, there are both standard Krylov estimators and randomized methods well-suited to this task. However, we establish barriers to solving this problem in the TF setting. In particular, if one uses Krylov queries, we show that the Arnoldi projected norm approximations to the spectral norm can exhibit nearly arbitrary convergence. We also prove a non-identifiability result for TF spectral and Frobenius norm estimation.

4.1 Krylov approximation of the spectral norm↩︎

We begin with the Krylov setting. The results in this section should be viewed in the spirit of the “any convergence curve is possible” literature for GMRES and Arnoldi, beginning with the classical result of Greenbaum, Pták, and Strakoš for GMRES convergence curves [11] and followed by related pathology results for restarted GMRES, Arnoldi Ritz values, and harmonic Ritz values [16]–[19]. Those works show that, under surprisingly weak constraints, Krylov output can be made to follow nearly arbitrary prescribed behavior. Our result is analogous in flavor but different in target: instead of residual curves or Ritz values, we consider the Arnoldi projected norms used as TF proxies for the spectral norm.

Let \(b\) of unit norm be given; we use this \(b\) as the first column \(v_1\) of \(V_k\). Recall that the Krylov space of order \(k\) is given by \(\mathcal{V}_k = \mathcal{K}_k(A,v_1) := {\rm span}\{v_1, Av_1, \dots, A^{k-1}v_1\}\). After \(k\) steps, the Krylov relation is \[A V_k = V_{k+1} \, H_{k+1,k}.\] Here, \(V_k\) and \(V_{k+1}\) are matrices with orthonormal columns, spanning \(\mathcal{V}_k\) and \(\mathcal{V}_{k+1}\) respectively, and \(H_{k+1,k}\) is an upper Hessenberg matrix. In a TF setting, a natural Arnoldi approximation to the spectral norm is the norm of the restriction of the operator to \(\mathcal{V}_k\), \[\|A_{\vert_{\mathcal{V}_k}}\| = \|AV_k\| = \|V_{k+1} \, H_{k+1,k}\| = \|H_{k+1,k}\|.\] (This quantity is at least as large as \(\|H_{k,k}\| = \|V_k^\top \! AV_k\|\).) The quantity \(\|H_{k+1,k}\|\) is the natural object produced by Arnoldi from forward matvecs alone. The point of this section is that, just as GMRES and Arnoldi can exhibit highly noninformative prescribed behavior in the classical pathology literature, these projected norms can also behave in an essentially arbitrary way. In particular, forward-only Krylov information need not provide a reliable guide to \(\|A\|\).

For the two results in this section, we use a technique that has also been used by [11], in the first “any convergence curve is possible” paper: a companion matrix which is also upper Hessenberg. The benefit of Hessenberg structure is that \(\mathcal{K}_k(A, e_1) = {\rm span}\{e_1, \dotsc, e_k\}\), while the companion matrix form makes it easier to prescribe the behavior of a specified desired quantity: an appropriate largest singular value (matrix two-norm) in our case. For consistency with the rest of the paper, our results are for matrices with real valued, however the identical proof also applies to complex matrices.

This proposition shows that the Arnoldi projected norms alone impose essentially no meaningful monotonic relation to the true spectral norm beyond the trivial bound \(\|H_{k+1,k}\|\le \|A\|\). In particular, a forward-only Krylov method may observe an arbitrarily prescribed nondecreasing history before the final norm is revealed.

The point of the next result is that the pathology from 2 is not merely a consequence of unconstrained spectral data. Even after fixing the eigenvalues, the projected norms produced by Arnoldi can still be forced to follow an arbitrary nondecreasing curve.

While the prior results exhibit arbitrarily bad convergence of a TF Krylov method, we investigate the typical behavior in the next two examples.

We now consider a family of finite difference matrices arising from PDEs.

4.2 Non-identifiability of a matrix norm↩︎

To estimate any matrix norm using matrix-vector products, we establish a lower bound of \(q \geq n\) matrix-vector products to ensure there is a unique solution to the problem. The gist of the result is that if one only has access to \(n-1\) matrix-vector products with the matrix \(A\), one can “hide” an arbitrarily large singular vector in a direction orthogonal to the \(n-1\) inputs. As a result, one can construct a matrix \(B\) satisfying the same \(n-1\) matrix-vector products as \(A\), but with an arbitrarily different norm. In particular, this result applies to the spectral, Frobenius, and Schatten norms.

This result highlights an important qualitative difference between query complexity lower bounds and non-identifiability results. While identifiability of the Frobenius norm requires \(n\) queries, one can use \(\mathcal{O}(\varepsilon^{-2})\) randomized, one-sided queries \(\|Ag_i\|^2\), where \(g_i\)’s entries are distributed as i.i.d. \(\mathcal{N}(0, 1)\), and apply Hutchinson’s estimator [54] to \(A^\top \! A\) to approximate \(\|A\|_F\) within relative error \(\varepsilon\). Because the identifiability result relies on an adversarial construction that depends on the transcript, a randomized algorithm can avoid this adversary with high probability.

5 Low-rank approximation↩︎

We now turn to low-rank approximation, a problem that is well-studied in the matrix-free context [55]–[59]. In particular, lower bounds on the two-sided query complexity of low-rank approximation have been thoroughly investigated [55], [56], [60]. In this section, we derive a simple lower bound on the query complexity of this task in the TF setting.

We see that low-rank approximation is one of the clearest examples in which TF access remains informative, but only at a prohibitive query cost. The key observation is that on exactly rank-\(k\) inputs, any \((1+\varepsilon)\)-approximation must recover the matrix exactly. This allows us to reduce approximate low-rank approximation to exact reconstruction on a carefully chosen hard distribution. We note that a non-identifiability result for this problem where the matrix is exactly rank-\(k\) is proven in [4]. Additionally, in the exactly rank-\(k\) case, an analysis of the space of possible row spaces or matrices that can be recovered from only forward matvecs is given in [61] and [2], respectively.

5.1 Query complexity lower bound↩︎

We prove a query complexity lower bound for the more general problem of finding a low-rank approximation that is within a constant factor \(\gamma > 0\) of optimal. Because this problem is strictly easier than relative approximation within a prescribed error \(\varepsilon> 0\), the lower bound also extends to relative approximation.

The theorem is sharp in a qualitative sense: in the TF model, near-optimal low-rank approximation is not a genuinely matrix-free task. On exact rank-\(k\) inputs, approximate recovery collapses to exact recovery, which requires full-dimensional information.

5.2 Low-rank approximation via column subset selection↩︎

The problem of column subset selection is well-studied in numerical linear algebra, particularly as it relates to matrix factorizations like interpolative/CUR decompositions [62]–[65], adaptive cross approximation [66], and the discrete empirical interpolation method (DEIM) [67]. We now discuss the problem of low-rank approximation via column subset selection, also known as finding a \(CC^\dagger\) or \(CX\) interpolative decomposition. The goal is to find an optimal or near-optimal set of \(k\) columns of \(A \in \mathbb{R}^{m \times n}\), stored as the columns of \(C \in \mathbb{R}^{m \times k}\) such that the low-rank approximation given by the projection \(CC^\dagger A\) satisfies \[\label{eq:css95lra} \|A - CC^\dagger A \|_F \le \gamma_k \cdot \|A - A_k \|_F.\tag{3}\]

Because one is restricted to using \(A\)’s columns for low-rank approximation, the theoretical best possible constant is \(\gamma_k \geq \sqrt{k+1}\) [68].

In the matrix-free setting of column subset selection, one may sample columns or rows by querying elementary basis vectors. Recently, both deterministic and randomized algorithms using matrix-vector products have been developed for matrix-free column subset selection with the guarantee 3  [69], [70]. However, one cannot form \(CC^\dagger A\) efficiently without transpose access, as one must apply \(C^\dagger\) to every column of \(A\). Hence, a key assumption in the randomized matvec algorithm in [69] is that in addition to forward matvec access with \(A\), one also has a good approximation to \(A\)’s row space, i.e., a matrix \(V \in \mathbb{R}^{n \times k}\) such that \(A\approx A VV^\top\). Prior results on TF row-space approximation already show that one cannot achieve this guarantee using only forward queries unless one performs \(n\) matvecs [2], [61].

If the task is just to find a good set of \(k\) columns for low-rank approximation without explicitly forming the low-rank approximation \(CC^\dagger A\), we can prove a non-identifiability result by reducing to the case where \(A\) is exactly rank-\(k\):

We explain in the following remark that the assumptions about \(X\) in 6 are necessary for non-identifiability to hold.

We have shown that, excluding two edge cases, low-rank approximation via column subset selection also fundamentally relies on transpose information, in the sense that the set of columns spanning the range of the matrix is not identifiable from fewer than \(n\) forward queries. In the following section, we address the related, more challenging problem of finding a maximum volume submatrix.

6 Local maximum volume submatrix↩︎

In this section, we consider identifiability for finding the local maximum volume submatrix in a TF manner. To describe the maxvol problem, for \(S\subset\{1,\dots,n\}\) with \(|S|=k\), write \[A_S = A_{:,S}, \qquad \operatorname{vol}(A_S) = \textstyle \prod_{\ell=1}^k \sigma_\ell(A_S) = \det(A_S^\top A_S)^{1/2}.\] The global maxvol problem is to find \[\max_{S \subset \{1, \dotsc, n\}, \, |S| = k} {\rm vol}(A_S).\] This problem is closely related to the problem of low-rank approximation via column subset selection described in the previous section. There are many contexts where one wishes to find a representative or dominant submatrix [71] and a variety of techniques using volume sampling/determinantal point processes [72] and leverage scores. While the global maxvol problem is NP-complete [73], it has been recently shown that the problem of finding the near-local maxvol of a matrix \(A\) is theoretically necessary for Gaussian elimination or the QR algorithm to be rank-revealing [74]. Moreover, this work proposes a simple greedy algorithm for this problem. The local maxvol solution also results in a choice of column indices that yield a near-optimal low-rank approximation, thus solving the column subset selection problem. However, an important distinction is that while column subset selection can achieve good low-rank approximation error in terms of the trailing singular values (as proven for, e.g., the methods of [70] and [69]), it does not necessarily guarantee a near-local maximum volume pivot or good singular value estimates.

The local maxvol problem is defined as follows. For a size-\(k\) set \(S\), define its one-swap neighborhood by \[\mathcal{N}_1(S) = \bigl\{(S\,\setminus\{i\})\cup\{j\}: i\in S,\;j\notin S\bigr\}.\] For \(\gamma\ge 1\), we say that \(S\) is a \(\gamma\)-local maximum-volume subset of \(A\) if \[\label{eq:gamma-local-maxvol} \operatorname{vol}(A_T) \le \gamma \cdot \operatorname{vol}(A_S) \qquad \text{for every } T\in \mathcal{N}_1(S).\tag{4}\] The case \(\gamma=1\) is ordinary one-swap local maxvol.

6.1 Identifiability of local maximum volume↩︎

We now consider identifiability for the local maxvol problem.

This proposition is a pure identifiability result: the issue is not computational efficiency, but the fact that the forward transcript does not determine the correct subset. Even exact optimization is impossible when distinct matrices agree on all queried products but induce different optimal column choices.

7 Schatten-\(p\) norm estimation↩︎

We next consider Schatten-norm estimation. For \(p \ge 1\), let \[\|A\|_{S_p} = \big( \textstyle \sum_i \sigma_i^p(A) \big)^{1/p}\] denote the Schatten-\(p\) norm. For even \(p\), we have \(\|A\|_{S_{p}}^{p} = \operatorname{tr}((A^\top \! A)^{p/2})\). Important examples of this quantity include the nuclear norm (Schatten-1), Frobenius norm (Schatten-2), and spectral norm (\(p \to \infty\)).

Schatten-\(2p\) norm estimation from samples is well-studied; in the context of covariance matrices, it is used via the method of moments to approximate the eigenvalues from observed samples drawn from the covariance matrix’s corresponding distribution [8]. This algorithm produces an unbiased estimator for the Schatten-\(2p\) norm using only one-sided matrix-vector products, and the variance of this estimator can be further reduced under additional assumptions, such as rapid decay of singular values [75]. In the case of Frobenius norm estimation (\(2p=2\)), this algorithm reduces to Hutchinson’s estimator applied to \(A^\top A\).

The query complexity of Schatten-norm and spectral norm estimation is studied in [9]. They consider the bilinear sketching model, where there is a distribution over \(r \times n\) matrices \(S\) and \(n \times s\) matrices \(T\) so that one observes \(SAT\) to approximate \(\|A\|\) within a constant factor. For even \(p \geq 4\), they obtain an \(\varepsilon\)-approximation to \(\|A\|_{S^p}\) with \(rs = \mathcal{O}(n^{2 - 4/p})\), which is also shown to be optimal in \(n\) and \(p\) dependence with a matching lower bound. Given this constraint on \(rs\), if one wants to minimize the total number of queries \(r + s\), the solution is \(r = s = \mathcal{O}(n^{1 - 2/p})\). We note that this result can be converted into a TF bound as follows. A non-adaptive TF algorithm using \(q\) forward matvecs in this form requires \(S\) to be invertible, so \(r = n\) and \(s = q\). Thus, substituting \(rs = nq\), a TF algorithm using non-adaptive, randomized queries of the form \(AT\) requires \(q = \Omega(n^{1 - 4/p})\) queries for even \(p > 4\). It is worth noting that the number of forward queries needed for spectral norm estimation, which can be recovered from \(\| A \|_{S^p}\) as \(p \to \infty\), is \(n\); see a related discussion in [76]. In this section, we obtain a more general lower bound for possibly adaptive queries. However, we lose the descriptive dependence on \(p\).

We focus first on Frobenius norm estimation, where transpose access leads to the classical contrast between Hutchinson’s estimator and Hutch++. Hutchinson’s estimator uses only one-sided information, while Hutch++ exploits a low-rank correction that relies on both \(A\) and \(A^\top\). The following theorem shows that this difference is not merely algorithmic: the improved \(\varepsilon^{-1}\) dependence is impossible in the TF model.

This theorem separates two phenomena that are often conflated. Hutchinson’s method is already optimal in the one-sided model, whereas the improved complexity of Hutch++ for estimating \(\|A\|_F\) fundamentally relies on access to both \(A\) and \(A^\top\). In the following remark, we distinguish this result from theory surrounding improved variants of Hutch++, which also address the trade-off between matvecs with \(A\) and \(A^\top\).

The following corollary extends 9 to general Schatten norms. Intuitively, as \(p \to \infty\), estimating \(\|A\|_{S_p}\) becomes harder because one requires more information about \(A\)’s spectrum. However, the following lower bound is independent of \(p\), and therefore very loose as \(p\) grows.

The corollary shows that, for lower bounds, the Frobenius norm already captures the difficulty of estimating every Schatten norm in the TF model. At least at the level of worst-case query complexity, the rank-one case is sufficient to transfer the obstruction to the whole Schatten family.

8 Least squares problems↩︎

We conclude the main body with a concrete non-identifiability result for overdetermined least squares problems; we note that a preliminary version of this result was included in the first author’s dissertation [7]. We consider the least squares problem \(Qx = b\), where \(b\) is known and \(Q\) is only accessed via forward matrix-vector products. We note this setting arises because \(A\) may be too large to store, but one can still store the right-hand side \(b\). Even in the simplest orthonormal setting, where the solution is just \(Q^\top b\), we show that forward matvecs alone do not determine the minimizer. This result complements recent lower bounds on the query complexity of approximately solving linear systems [10]. Let \(\kappa(B) = \sigma_{\max}(B) \, / \, \sigma_{\min}(B)\) be the condition number of \(B\), where \(\sigma_{\max}\) and \(\sigma_{\min}\) are the largest and smallest singular values of \(B\), respectively.

In the following remark, we explain that the assumption that \(x \in \text{col}(X)\) in 11 is necessary to for non-identifiability.

Thus, even in the benign orthonormal setting, forward access to the design matrix does not determine the least squares minimizer. The obstruction again comes from an unseen direction: one can perturb the problem without changing any queried products, while moving the solution and keeping the matrix well conditioned. The distance between the solution \(x\) and the query space \(\text{col}(X)\) determines the identifiability of the least squares solution.

9 Conclusions↩︎

We have shown that one-sided access to a matrix is often substantially weaker than standard matrix-free access with both \(A\) and \(A^\top\) available. On the pathology side, Arnoldi projected norms can behave arbitrarily, so forward-only Krylov information does not by itself provide a trustworthy proxy for the spectral norm. On the complexity side, near-optimal low-rank approximation requires \(n\) forward matvecs, while Frobenius norm estimation requires \(\Omega(\varepsilon^{-2})\) queries when \(n\) is sufficiently large. On the identifiability side, solutions to the local maximum volume submatrix problem and orthonormal least squares are not determined by forward transcripts alone. Taken together, these results suggest that transpose access is not merely a technical convenience, but a structural source of information.

This does not necessarily mean that the TF setting is hopeless. Rather, it suggests that when TF algorithms succeed, there must be some additional structure in the matrix at hand. For example, while adjoint-free operator learning may be impossible or highly inefficient for general operators, solution operators of elliptic PDEs provide indirect information about the adjoint through forward queries [2]. In other cases, one may have a priori knowledge about the operator’s row space [69], or noisy interactions with the transpose operator in the “unmatched transpose” setting. Moreover, our experiments in 1 and  2 demonstrate that for matrices exhibiting singular value decay or arising from PDEs, TF methods can even outperform ones which use the transpose. Thus, a natural continuation of the present work is to move beyond worst-case analyses and understand the performance of TF methods for structured matrices that commonly arise in practice.

Many other natural directions remain open for exploration. Are there stronger lower bounds for Schatten-\(p\) norm estimation beyond the rank-one reduction used in 1? In particular, can one extend the \(p\)-dependent lower bounds in [9] to adaptive sketches? What is the correct TF complexity of singular-vector or nullspace estimation? Moreover, to what extent can analogous impossibility results be established for regularized least squares problems or other structured inverse problems? We hope that the present work provides a starting point for these deeper questions in TF linear algebra.

Acknowledgments↩︎

D. H. thanks the Simons Institute for the Theory of Computing for hosting her as a research fellow in the Fall 2025 program on Complexity and Linear Algebra. She also wishes to acknowledge many useful conversations with Cameron and Chris Musco on lower bound techniques and Pinsker’s inequality. A. T. has been supported by the Defense Advanced Research Projects Agency (DARPA) through The Right Space (TRS) Disruption Opportunity (DARPA-PA-24-04-07). The views, findings, and conclusions expressed in this paper are those of the authors and do not necessarily reflect the official policy or position of DARPA, the U.S. Department of Defense, or the U.S. Government. A. T. has also been supported by the National Science Foundation CAREER grant DMS-2045646.

References↩︎

---

1. Courant Institute, New York University, New York, NY 10012, USA, diana.halikias@nyu.edu.↩︎

2. Department of Mathematics and Computer Science, TU Eindhoven, The Netherlands, m.e.hochstenbach@tue.nl.↩︎

3. Mathematics Department, Cornell University, Ithaca, NY 14853-4201, USA, townsend@cornell.edu.↩︎

4. Date 2026-06-12↩︎

5. We consider algorithms that use \(\mathcal{O}(n)\) forward queries, where \(n\) is the number of columns in the matrix, to be inefficient. This contrasts with the limited memory setting, which requires \(\Omega(n)\) matvecs [6].↩︎

6. In fact, [14] was the most cited math paper of the 1990s.↩︎

7. We gratefully thank Chris and Cameron Musco for their help on this topic.↩︎