Skip to main content

Questions tagged [matrix]

For questions about using and representing matrices on a computer in order to solve computational problems. Should generally also include a tag about the specific property/problem you are solving (e.g. [tag:linear-algebra], [tag:eigenvalues], [tag:inverse].

Filter by
Sorted by
Tagged with
10 votes
1 answer
319 views

Wanted: sequences of linear systems for recycling Krylov solver analysis

In the solution of sequences of linear systems $$A_ix_i=b_i\quad\text{for}\quad i=1,2,\dots$$ with Krylov subspace methods, data can be recycled from already solved linear systems in order to speed up ...
2 votes
0 answers
56 views

What is fastest method for finding the minimum and maximum eigenvalues of a (possibly very large) symmetric matrix?

What is the best way to find the extreme eigenvalues - in order to find the spectral radius - of a general real dense symmetric matrix? Looking at similar questions e.g.: What's the most efficient ...
2 votes
1 answer
60 views

Is there any matlab built-in function or libraries to calculate $\frac{d(\ln A)}{dA}$?

we can first conduct spectral decomposition of an positive definite isotropic tensor $A$ and then we can define $\ln(A)$, then we can define the frechet derivative of it, but how to calculate this in ...
5 votes
1 answer
115 views

Optimized Lanczos method for finding eigenvalues of $A \otimes B$

Recently my supervisor told me about an efficient way to calculate eigenvalues and eigenvectors of matrix $A \otimes B$ with $a_{1} \times a_{2}$ as dimensions of $A$ and $b_{1} \times b_{2}$ is of $B$...
0 votes
0 answers
31 views

Lumped (diagonal) vs. consistent (non-diagonal, symmetric) mass matrix in Nastran

I've been tinkering with DMAP to explore the procedure followed by Nastran when solving a complex modes analysis. I've reached a passage I cannot understand: at some point Nastran formulated what it ...
0 votes
1 answer
74 views

What do diagonal (DOF-to-self) terms of stiffness matrix physically mean?

I am used to interpreting each entry of a solid mechanic system's stiffness matrix as a 1D (linear or angular) spring joining one DOF (column index) to another (row index). But this interpretation ...
1 vote
0 answers
57 views

Orthogonal Transformation of Hessenberg Matrices

$H\in\mathbb{R}^{n\times n}$ is an upper Hessenberg matrix. Suppose $\lambda$ is an eigenvalue of $H$ and $x$ is an eigenvector w.r.t. $\lambda$. Is there any fast algorithm that can find an ...
4 votes
2 answers
249 views

Reordering eigenvalues in Schur factorization - MATLAB ordschur and LAPACK dtrsen not producing the same results

Disclaimer: I previously posted this on SO, but though it would be more relevant for scicomp. The original post has been deleted. I have been trying to recreate the functionality provided by MATLABs <...
2 votes
2 answers
170 views

Confusion about matrix differentiation in a nonlinear matrix equation

I am trying to solve a matrix equation in the following discrete form: $$ \frac{K^{n+1}-K^n}{\Delta t} = [(K^{n+1} (V^{n})^T).^3 - K^{n+1} (V^{n})^T]V^n. $$ where $K^{n+1} \in \mathbb{R}^{m \times r}, ...
0 votes
1 answer
68 views

Weird runtime behavior of `scipy.linalg.solve_triangular` and `trtrs`

I want to understand the time complexity of scipy.linalg.solve_triangular, which calls trtrs from LAPACK under the hood, so I ...
1 vote
0 answers
107 views

How to implement boundary conditions for the Thomas algorithm

For my variable $U(t,x)$, I have implemented the thomas algorithm with $U_j^i$: $$ a(x)U_{j-1}^{i+1}+ b(x)U_j^{i+1} + c(x)U_{j+1}^{i+1} = d(x)U_j^{i} $$ Then $\textbf{A}$ is a tridiagonal vector with ...
1 vote
1 answer
189 views

How do BLAS libraries implement support for transposed matrices?

I'm trying to understand how BLAS libraries implement fast GEMM with support for transposed matrices. Say, I'm only operating on square matrixes (with dimensions n ...
1 vote
0 answers
124 views

How to vectorise numerical differentiation

I have a 2-D matrix with 2 spatial coordinates and I want to be able to vectorise the process of numerically differentiating with respect to its 2 coordinates, rather than just looping along the rows ...
5 votes
2 answers
584 views

$\mathbf{UDU}^\top$ decomposition routines in LAPACK/Eigen?

I would like to compute the decomposition of a real symmetric positive definite matrix $\mathbf{A} = \mathbf{UDU}^\top$. LINPACK seems to have it as DSIFA, but I ...
11 votes
1 answer
342 views

Is it possible to express an arbitrary tensor contraction in terms of BLAS routines?

I noticed that libraries like numpy and pytorch are able to perform arbitrary tensor contractions at speeds similar to comparably sized matrix multiplications. This leads me to believe that underneath ...
3 votes
1 answer
92 views

Apply 3D Operator to Matrix and get new Matrix

Hopefully this question makes sense. I know I can formulate an operator for a vector as a matrix, then apply that matrix to my vector to get a new vector. For example, if I define a left shift ...
0 votes
0 answers
69 views

Singular Matrix Error in Incomplete LU Decomposition

I’m currently working on solving the following PDE: $$\begin{equation} -(\mu_x \frac{\partial^2 u}{\partial x^2} + \mu_y \frac{\partial^2 u}{\partial y^2}) = f(x, y)\end{equation}$$ Where a right hand ...
5 votes
2 answers
2k views

What algorithm(s) do numpy and scipy use to calculate matrix inverses?

I am solving differential equations that require inverting dense square matrices, and I wanted to know what algorithm(s) do numpy and scipy use to calculate matrix inverses?
3 votes
1 answer
149 views

Overlap matrix and its inverse matrix

Now, we consider a non-orthonormal basis: $$\mathcal{S}_N=\{|\alpha\rangle,a^\dagger|\alpha\rangle,a^{\dagger 2}|\alpha\rangle,\ldots,a^{\dagger N}|\alpha\rangle\},$$ where $|\alpha\rangle$ is the ...
0 votes
1 answer
108 views

Finding the Vector $v$ for a Given Householder Matrix Transformation of Non-Collinear Vectors $a$ and $b$

Consider a vector $v$ in $\mathbb{R}^{n\times1}$. The Householder matrix is defined as follows: $$H(v)=I-\dfrac{2vv^T}{v^Tv}.$$ It can be demonstrated that $H(v)$ is symmetric and orthonormal. The ...
9 votes
1 answer
1k views

Why is it that SVD routines on nearly square matrices run significantly faster than if the matrix was highly non square?

In Python / Matlab, if you run a routine for SVD on a significantly non-square matrix, X, such as X.shape = (2,15000) you will ...
2 votes
0 answers
67 views

Efficient Algorithm for LU-Factorization of Modified Matrix with Last Column Alteration If We Have Its Not-Modified LU-Factorization

Suppose that we have a $n\times n$ matrix $A$. We have its LU-factorization as $A=LU$ (or $PA=LU$ that $P$ is a permutation matrix). Now assume we change the last column of matrix $A$ and denote the ...
3 votes
1 answer
173 views

How can we calculate mixed derivatives numerically using the Chebyshev derivative matrix?

Using the Chebyshev derivative matrix $D$, we can numerically approximate the first and second derivative of a function by doing matrix multiplication: $${df(x) \over dx} = Df(x) \tag 1$$ $${d^2f(x) \...
2 votes
0 answers
73 views

Complex matrix logarithm discontinuity by solving inverse Fourier integral by alternative method to FFT

NOTE: This code is a piece of code I am using for a master's thesis, so I do not expect someone to do the work for me, but I gladly accept suggestions of any kind. However, I am trying to get the ...
16 votes
2 answers
14k views

In FEM, why is the stiffness matrix positive definite?

In FEM classes, it's usually taken for granted that the stiffness matrix is positive definite, but I just can't understand why. Could anyone give some explanation? For instance, we can consider the ...
1 vote
0 answers
79 views

numerical calculation of haldane model arm chair edge states

hello I am trying to numerical simulate the band structure of the one-dimensional periodic arm chair edge states, I use the pybinding model to construct and ...
1 vote
0 answers
36 views

Sums of sparse matrices modulo 2

I have a set of sparse matrices modulo 2 (i.e systems of equations in the format "a XOR b XOR ... XOR c = 0 or 1" compressed by me into reduced row echelon form). This is a rather large set (...
2 votes
1 answer
134 views

Tools to compare two matrices with same dimensions

Context: I have two 3D non-random matrices that have the same dimensions. These matrices represent satellite images with 1 band, so their values are strictly positive. They both present areas that ...
2 votes
1 answer
178 views

Measuring the extent to which two sets of vectors span the same space

I have a set of measurements $y_i$, $1 \leq i \leq N$, and I want to model these measurements with a linear model. I have two possible models I can use, $$ y \approx A c $$ and $$ y \approx B d $$ ...
1 vote
1 answer
128 views

Calculating camera calibration matrix with Scilab

I'm not entirely sure whether this question belongs here or in DSP but I think this is the proper site. I'm following these videos (first video, second video) to calibrate a camera for photogrammetry ...
1 vote
1 answer
95 views

Automatic differentiation (AD) of a loss function which maps unitary matrix onto number

Is it possible to estimate whether automatic differentiation (AD) techniques could enable a more efficient way to repeatedly compute the derivative $\delta L / \delta u^*_{ij}$ of a specific loss ...
1 vote
1 answer
85 views

Reshaping a matrix and rearranging the elements lexicographically into a vector

Let us say that I have a $3 \times 3$ matrix $\bf X$, that has to be reshaped into a vector and rearranged as follows: $$ {\bf v} ({\bf X}) = \begin{bmatrix} x_{22} & (x_{23}+x_{32}) & (x_{21}+...
5 votes
1 answer
544 views

Converting distance matrix back into original data

Suppose that we have $N$ points, and a distance matrix $D \in \mathbb{R}^{N \times N}$ describing the Euclidean distance among those points. For now, assume that we do not necessarily know how many ...
2 votes
0 answers
73 views

When is Lanczos tridiagonalization accurate?

Suppose that we are given a random, symmetric matrix $A$, and a random vector $q$. For concreteness, assume the dimensions of $q$ and $A$ are both $1,000$. I would like to use the Lanczos algorithm to ...
1 vote
1 answer
191 views

QR algorithm for eigenvalues and eigenvectors of large symmetric matrices

I am trying to write a QR algorithm in Python for eigenvectors and eigenvalues finding for large symmetric matrices, My initial thought was to use Householder transformation with a Wilkinson shift on ...
2 votes
1 answer
483 views

In Octave, how do I specify that the solution to a matrix equation should be over integers?

In Octave, how do I specify that the solution to a matrix equation should be over integers? I.e., Given matrix $A$, vectors $x$ and $b$; $Ax=b$. Find vector $x=A^{-1}b$ such that all its entries are ...
2 votes
1 answer
224 views

Parallelize pseudo inverse of a matrix using Lapacke

I am currently using the protocol described in https://stackoverflow.com/questions/55599950/computation-of-pseidoinverse-with-svd-in-c-using-blas-and-lapacke to compute the pseudo inverse of a matrix. ...
4 votes
5 answers
654 views

Matrix derivative

I am looking to compute the derivative of the following expression: $$\frac{\partial}{\partial X}\mathrm{tr}\left[A\exp(X)\right]$$ where $A$ is both a symmetrical and positive-definite matrix and $X$ ...
3 votes
1 answer
123 views

Scipy solve_ivp sensitivity to random phase shifts

I am trying to solve a coupled system of ODE's using the solve_ivp function from scipy. The general form of the equation is given via $$\dot{y}(t) = M(t)y(t).$$ The time dependence of matrix is ...
19 votes
3 answers
3k views

Why do we usually not want the eigenvalues of non-symmetric matrices?

I came across this line in a class note I am reading where it discusses finding eigenvalues of matrices. In reality we don't go all the way with Arnoldi. We stop at a decent value of 𝑘. Then the 𝑘 ...
4 votes
5 answers
326 views

Determining the algorithmic complexity

A few of the iterative matrix algorithms (CG,GMRES etc.) I have authored are acting rather funny. They converge to the right answers but take abnormally long time to run. I am in the process of ...
0 votes
1 answer
203 views

How to efficiently transpose distributed matrix in Scalapack?

I have a distributed matrix in block cyclic layout. Is there an efficient way to out/in place transpose a distributed matrix with scalapack? Context: I am trying to diagonalize the transpose of a ...
3 votes
1 answer
141 views

Rewriting quadratically-constrained optimization problem as a semidefinite program

Suppose $A,H$ are positive definite matrices and $\alpha,t$ are scalars. Is there a way to massage the following problem into a form suitable for a specialized solver? $$\begin{array}{ll} \underset{\...
1 vote
1 answer
485 views

2-norm and infinty norm of a system in controls

How to compute 2-norm or infinity norm of following system? i am confused whether to calculate using simple matrix theory "where it don't regard for s domain" or H2 and H-infinty norm. ...
1 vote
0 answers
135 views

Does cblas_dgemm mutate my input matrices?

I have written a matrix class Matrix<T> for which I have implemented a wrapper function for cblas_dgemm. ...
0 votes
3 answers
294 views

How do we compute a matrix-vector product without using the standard matrix multiplication?

I am a bit confused right now. I am taking a class on numerical linear algebra and as I have understood sometimes one doesn't compute the Matrix-Vector product $Av$ normally but uses other techniques ...
1 vote
1 answer
142 views

QR decomposition with only diagonal elements changing

I have to compute the QR decomposition of a matrix A repeatedly, each time with ONLY diagonal elements changing. Is there an efficient way to accomplish this ...
3 votes
1 answer
273 views

BiCGSTAB convergence

So I need a fast converging solver for SysLinEq as a subroutine in fortran, decided to test BiCGStab in Matlab. Thank God I decided to test it out on first before implementing in Fortran as a ...
2 votes
1 answer
516 views

Matrix regularisation for ill-conditioned problems

I have read that matrix regularisation can improve the stability of LU or Cholesky decomposition of ill conditioned problems. The idea is to add a value to the diagonals of a matrix: $B=A+cI$ In the ...
2 votes
0 answers
57 views

Is there a permutation used in sparse QR factorizations that better locates small elements on the diagonal of R?

Is there a permutation used in sparse QR factorizations that better locates small elements on the diagonal of R to the end of the diagonal? As an example, consider the following snippet of MATLAB ...

1
2 3 4 5
13