Questions tagged [linear-algebra]

Questions on the algorithmic/computational aspects of linear algebra, including the solution of linear systems, least squares problems, eigenproblems, and other such matters.

Filter by
Sorted by
Tagged with
0
votes
1answer
47 views

Fast iterative approximate order-oblivious Orthogonalization algorithm?

I have set of N m-dimensional vectors $\{\phi_i\}$ which gradually loose mutual orthogonality in an algorithm. => I have to re-orthogonalize them every few iterations. But if I do e.g. Gram–Schmidt ...
3
votes
1answer
57 views

Why are the round-off errors when solving the linear system $Ax = b$ of order $\varepsilon_\text{mach} x_j$?

I was reading a paper on arXiv where, in Section 2.4, the authors are discussing the error that arises in the solution of a linear system $$Ax = b,$$ or, to match up better with the paper, $$\Phi \...
4
votes
0answers
45 views

Is there a library that allows einstein summation on dense, sparse, and LinearOperator type tensors

Numpy's einsum only works with dense tensors. Is there an alternative that also works with sparse tensors and linear operators? For example, I might have a ...
3
votes
1answer
95 views

Matrix condition number and reordering

Does the condition number change when a matrix is reordered by e.g. Cuthill Mckee or some other method?
4
votes
0answers
68 views

Is there any catch on using `zgemm3m` vs regular `zgemm`?

I've just (to my embarrassment) encountered a BLAS-like extension of a matrix-matrix product subroutine gemm in Intel MKL: gemm3m...
1
vote
1answer
61 views

How to implement this even-odd matrix decomposition efficiently?

Note: This question has also been asked on stackoverflow - see https://stackoverflow.com/questions/57197910/how-to-implement-this-even-odd-matrix-decomposition-efficiently?noredirect=1#...
1
vote
0answers
40 views

Does Boost provide a template implementation of the wedge product?

Does the boost C++ library implement the computation of the wedge product? The wedge product is mentioned here, but it is not very clear (to me at least) whether there is a template implementation of ...
1
vote
0answers
73 views

What's the more efficient way to solve this matrix equation?

This is intended to be a more generic question not about a specific system. Given a hermitian matrix $H(x_1,\dots,x_n)$ depending non-linearly on some real parameters $x_1,\dots,x_n$. We want these to ...
2
votes
1answer
36 views

How to pass matrices to parallel workers quickly in matlab?

I am trying to solve many different linear systems in parallel in matlab. The problem is, each linear system has entirely different parts and are fairly large, so passing the information to each of ...
1
vote
1answer
73 views

Sparse matrix inversion

I have the impedance matrix $Y$, formulated from an electrical network by augmented nodal analysis. The matrix $Y$ is shown as an image to illustrate its feature visually, where all the white blocks ...
2
votes
1answer
34 views

numpy.outer without flatten

$x$ is an $N \times M$ matrix. $y$ is a $1 \times L$ vector. I want to return "outer product" between $x$ and $y$, let's call it $z$. z[n,m,l] = x[n,m] * y[l] ...
0
votes
1answer
47 views

what does “D = diag(W.1)” means?

, what does “D = diag(W.1)” means?on page #2, just below equation (6) PFA screenshot and here is the link of the paper - original paper
4
votes
1answer
63 views

Fast algorithm for computing cofactor matrix

I wonder if there is a fast algorithm, say ($\mathcal O(n^3)$) for computing the cofactor matrix (or conjugate matrix) of an $N\times N$ square matrix. And yes, one could first compute its determinant ...
0
votes
0answers
46 views

Null space for smoothed aggregation algebraic multigrid

I do not really get the point of null space usage for creating the prolongation operator for smoothed aggregation algebraic multigrid. I know what the null space is per definition and I know that the ...
2
votes
0answers
66 views

Regularized least squares with QR factorization

Consider the regularized least squares problem $$ \min_x || b - A x ||^2 + \lambda^2 ||x||^2 $$ which is equivalent to $$ \min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{A \\ \lambda I} x \right|\...
4
votes
0answers
64 views

Complex differentiation of linear solvers

I have a linear system $$Ax=b$$ which I'm solving approximately, and I need to take the frechet derivative of x with respect to z. Were I solving the problem exactly (either analytically or to machine ...
3
votes
1answer
83 views

Optimality of block-Jacobi preconditioner

For a dense $N \times N$ matrix $A$, is the block-Jacobi preconditioner comprising the inverse of the diagonal blocks of $A$ the optimal block-diagonal preconditioner? Could there exist another matrix ...
1
vote
0answers
46 views

Kinetic preconditioning

Publication arXiv:0804.2583 describes a method for doing self-consistent iteration without having to diagonalize the Hamiltonian operator at every step. IX. PRECONDITIONING As already ...
0
votes
0answers
26 views

A preconditioner for self-consistent iteration

I tried to derive a preconditioner for self-consistent iteration similar to section IX in arXiv:0804.2583. For simplicity, consider here only one orbital (one or two electrons) systems. Suppose that ...
0
votes
0answers
68 views

Algebraic multigrid for coupled equations

As far as I understand is algebraic multigrid(AMG) a method that was intentionally developed to solve linear systems where every grid point or node has a single DOF. When AMG should now be used for ...
3
votes
1answer
63 views

Whitening transformation does NOT return a unit covariance matrix

For this question, I am using the following Wiki definition of Matrix whitening: Suppose $X$ is a random (column) vector with non-singular covariance matrix $\Sigma$ and mean 0. Then the ...
2
votes
1answer
139 views

In iterative methods, are matrix decompositions considered useful for implementation?

When we study an iterative method from textbooks, for example, see the Gauss-Seidel Method, the given matrix is decomposed with suitable splittings. In the example, $A = L+U$. So we can proceed with ...
5
votes
2answers
349 views

Python: vectorizing a structured linear system solve

Overview I am looking for a way to solve a structured linear system in Python without using a for loop (preferably using vectorization, if possible). Background Consider the following linear system:...
2
votes
1answer
66 views

Re-using LU factorization within iterative (?) setup for a sum of two matrices

So, I would love to make at least some use of my preexisting data, no matter how small, and just out of ideas. Maybe I am just a prisoner of a Kahneman-like theatre-ticket paradox, and don't know ...
4
votes
1answer
153 views

Accurate way of getting the square root inverse of a positive definite symmetric matrix

What is the most accurate algorithm to get the square root inverse of a positive definite symmetric matrix? I am not looking as much for efficiency, though using quadruple precision computation is out ...
2
votes
2answers
175 views

Compiled c++ code runs much faster with double than float. Explanation?

I am still rather new on here and I hope question is suitable for this forum otherwise please help me migrate it to greener pastures. I am an electrical engineer specializing in applying mathematics ...
5
votes
1answer
79 views

CHOLMOD condition number estimate

The CHOLMOD library provides a CHOLMOD_rcond function that estimates the reciprocal condition number (in the one norm) of a symmetric positive definite matrix from ...
4
votes
2answers
70 views

How to efficiently invert $K \otimes M+I_T\otimes \Sigma$?

I'm looking for a way to efficiently invert $$K \otimes M+I_T\otimes \Sigma$$ where the inverses for $M,K$ exist. $I_T$ is the identity matrix of dimension $T$, and $\Sigma$ is a diagonal matrix, with ...
2
votes
1answer
36 views

Discretization with non-constant matrix containg entries form unknown vector

Consider a system of PDEs $$ \begin{cases} u_t = \nabla \cdot (D(u)\nabla u) + \frac{c}{K_U+c}u-ku\\ c_t = d_c\Delta c -\frac{\nu_U c}{K_U + c}u \end{cases} $$ with some boundary conditions. Here, $D(...
29
votes
16answers
9k views

Good examples of “two is easy, three is hard” in computational sciences

I recently encountered a formulation of the meta-phenomenon: "two is easy, three is hard" (phrased this way by Federico Poloni), which can be described, as follows: When a certain problem is ...
1
vote
0answers
37 views

Stability of SVD, Eigendecompositions, and pseudoinverse procedures in modern LAPACK routines

I have proposed an optimisation algorithm which I claim has improved upon the previous algorithm in a number of ways. One of these claims is that my proposed solution requires no explicit SVD and ...
2
votes
1answer
79 views

How to add extra constraints to a linear system for probabilities?

Background: I have an equation which looks like as follows: $W \times P = R$ $$\left[\begin{array} &{1}&{0}&{0}&-\frac{w_{1}}{w_{o1}} &\dots &{0} &-\frac{w_{1}}{w_{0} } \...
2
votes
1answer
61 views

Rank of Hadamard Product with Masked Matrix

I have a matrix $A\in\{0,1\}^{d\times n}$ and $rank(A)=d,d<n$, and another matrix $X\in \mathbb{R}^{d\times n}$, but I do not know the rank of $X$. What can we say about the rank of their Hadamard ...
2
votes
1answer
75 views

How to compute the determinant of Hessian of a multivariable function?

I have a function $F(\vec x)$ of many variables (let's say in the order of hundreds of thousands). I need to compute the determinant of the Hessian matrix at the point $x_0$. Is there a way to ...
0
votes
0answers
54 views

Bound for Expectation of Singular Value

In my case, $X_{\boldsymbol{\delta}}\in\mathbb{R}^{d\times M}$ is a function of Rademacher variables $\boldsymbol{\delta}\in\{1,-1\}^M$ with $\delta_i$ independent uniform random variables taking ...
3
votes
2answers
67 views

Convexity of Sum of $k$-smallest Eigenvalue

If I have a real positive definite matrix $A\in\mathbb{R}^{n\times n}$, and denote its eigenvalues as $\lambda_1\leq \lambda_2 \leq ... \leq \lambda_n $. Define the function as $f(A)=\sum_{i=1}^{k} \...
0
votes
1answer
46 views

An Upper Bound of $\left<H,X\right>_F$ with Constrainted Rank

We have $X=[\mathbf{x}_1,\mathbf{x}_2...,\mathbf{x}_n]\in\mathbb{R}^{d\times n}$, $H=[\mathbf{h}_1,\mathbf{h}_2...,\mathbf{h}_n] \in\mathbb{R}^{d\times n}$, and $d<n$. $H$ has rank $r\leq d$ and $X$...
2
votes
1answer
82 views

Which pseudo-inverse to compute when Inverse is not possible? (No linear solve)

Let us assume that we have a function, $f(A)=\text{vec}(A^{-1})^\intercal B$, dependent on $A^{-1}$. However, due to some machine-precision limitations, the programming language I'm using cannot ...
5
votes
2answers
173 views

Computing any element of the null space of a singular matrix

Given a singular matrix $A$, what is the fastest method to find a single non-zero solution to $Ax=0$? Note that we are not looking for the whole kernel, we just want any non-zero vector in it. I ...
2
votes
1answer
109 views

A fast way to check if a Matrix is ill-conditioned, and turning it into well-conditioned

I'm running a simulation, and some linear solvers are returning a message of ill-conditioned matrix. Hence, I'm looking for a fast, easy to implement, method to detect if a matrix is ill-conditioned, ...
2
votes
1answer
22 views

Determine image of hypercube under linear map

Let $A$ be an $3\times N$ matrix (where $N$ is large) with nonnegative real entries. I'd like an algorithm for determining when a vector $v\in\Bbb R^3$ can be written as $Aw$ for some vector $w\in\Bbb ...
1
vote
0answers
93 views

Connection between piecewise linear basis functions and RELU activation function

ReLU activation is defined as follows $$\sigma(x)=\max(0, x).$$ Let's assume that I have deep network of 1 hidden layer, than output from my layer has form $$ f(x)= \sigma(Wx +b), $$ where matrix W ...
1
vote
1answer
68 views

Gaussian Elimination Using Fortran [closed]

I developed the code below for performing gaussian elimination in order to evaluate the determinant of a matrix: ...
1
vote
1answer
61 views

Product of rank one updates as a low rank update for quasi newton/BFGS

I'm trying to improve the speed of the following iteration to calculate $s_k$: $$B_k^{-1} = \Bigg( I + \frac{s_{k}s_{k-1}^T}{||s_{k-1}||^2}\Bigg)...\Bigg(I+ \frac{s_1s_0^T}{||s_0||^2}\Bigg) B_0^{-1}\\...
9
votes
3answers
331 views

fastest linear system solve for small square matrices (10x10)

I am very interested in optimizing the hell out of linear system solving for small matrices (10x10), sometimes called tiny matrices. Is there a ready solution for this? The matrix can be assumed ...
0
votes
1answer
56 views

Finding Matrix inverse with LU and repeted left division calls

Hello I am in a basic numerical methods class and our teacher has given us an algorithm which can compute the inverse of a matrix other than using MATLAB's built in library function. ...
5
votes
3answers
79 views

Finding exact rational solution to linear integer equations in Matlab

I have a linear system of equations $$Ax=b$$ where $A$ is an $N\times N$ matrix with integer values, and $b$ is a $N\times 1$ vector with integer values. Due to prior knowledge, I know that I am ...
6
votes
1answer
108 views

numerical solution of an under-determined linear equation in high dimensions

I need to solve a linear regression problem $$Ax=y$$ which is hugely underdetermined. I have around $10^6$ features but only $10^3$ equations. So $A$ is a $1,000\times 1,000,000$ matrix and $y$ a ...
2
votes
1answer
104 views

Projection onto the set of Orthogonal matrices

Let $M \in \mathbb{R}^{n \times n}$ and denote the set of Orthogonal matrices by \begin{equation} \mathcal{O}_{n} = \left\lbrace Q \in \mathbb{R}^{n \times n} \colon QQ^{T} = \mathbb{I}_{n} \right\...
2
votes
2answers
115 views

Asymptotic Complexity of Gaussian Elimination using Complete Pivoting

I would like to know the algorithm asymptotic complexity with Complete Pivoting. With partial pivoting, it is known to be $O(n^3)$. Is it the same for complete pivoting?