Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I have an application where I have to minimize a cost function of the form:
$J(\mathbf x) = \| \mathbf A \mathbf x - \mathbf b \|^2 \quad (1)$

By calculating the gradient, I derived that I have to solve the system of equations:
$\mathbf A^T \mathbf A \mathbf x = \mathbf A^T \mathbf b \quad (2)$

Now my question is, when can I solve the following system instead?
$\mathbf A \mathbf x = \mathbf b \quad (3)$

From my point of view, this depends on $\mathbf A$ to be invertible. In my application $\mathbf A$ is a square matrix of the form $\mathbf A = \mathbf I - \mathbf W$ where $\mathbf I$ is the identity matrix and $\mathbf W$ is a square matrix with zeros on the main diagonal and small values on a few secondary diagonals. The values of $\mathbf W$ are arbitrary but normalized, so that each row of $\mathbf W$ sums to 1. However, some lines of $\mathbf W$ can also be zero.

For example $\mathbf A$ may look like this:
$\mathbf A = \pmatrix{ 1 & -0.2 & 0 & 0 & -0.3 & -0.2 & -0.3 & 0 & 0 \cr -0.1 & 1 & -0.2 & 0 & 0 & -0.4 & -0.2 & -0.1 & 0 \cr 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \cr 0 & 0 & -0.2 & 1 & -0.1 & 0 & 0 & -0.6 & -0.1 \cr -0.2 & 0 & 0 & -0.2 & 1 & -0.5 & 0 & 0 & -0.1 \cr -0.4 & -0.3 & 0 & 0 & -0.1 & 1 & -0.2 & 0 & 0 \cr -0.1 & -0.1 & -0.1 & 0 & 0 & -0.6 & 1 & -0.1 & 0 \cr 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \cr 0 & 0 & -0.2 & -0.4 & -0.1 & 0 & 0 & -0.3 & 1 \cr }$

Have I argued correctly? Then how can I show that $\mathbf A$ is invertible? Or is there any other argument for solving (3) instead of (2)?

I tried to solve both systems (2) and (3) with MATLAB and the Intel MKL, but surprisingly only (3) gave me the expected result. I would expect that it also works with (2). Maybe a numerical problem?

share|improve this question
    
What are the dimension of your matrix and the vectors $x$ and $b$. –  Mhenni Benghorbal Jul 20 '12 at 11:44
    
It's actually an image processing problem, where the dimension $n$ equals the number of pixels. (e.g. $n=2073600$ for a high definition image). $\mathbf A$ is a $n$x$n$ matrix. $\mathbf x$ and $\mathbf b$ are $n$x1 column vectors. –  benni Jul 20 '12 at 13:19
1  
FYI, you might want to use conjugate gradient to solve $A^TA x = A^Tb$ rather than solve $Ax=b$ directly. In image processing and the like, you often want an iterative method that applies the operator without actually forming the matrix. It's not necessarily the fastest algorithm since you're roughly squaring the condition number, but in my experience it's easy to implement and get an answer quickly. –  p.s. Jul 20 '12 at 15:39
    
I also thought about doing that. However, since I don't get any useful result when solving $\mathbf A^T \mathbf A \mathbf x = \mathbf A^T \mathbf b$ using a direct solver, I haven't tried any iterative solvers yet. –  benni Jul 20 '12 at 16:17
    
Note: $A^TAx=A^Tb$ is more well behaved and perhaps easier to solve. It is symmetric and positive definite (or perhaps semi-definite). –  Gabriel Landi Jul 21 '12 at 1:21

4 Answers 4

Well if there exists an $x$, so that $Ax = b$, then of course the minimum of $J(x)$ will be $0$ (because $J(x) \geq 0 ~\forall x$). There exists a solution $x$, if and only if rank(A) = rank(A,b). So it doesnt depend on A being invertible, rather if b "adds" to the rank of A.

If there exists no such $x$, you have to solve $A^TAx = Ab$

share|improve this answer

Just some thoughts: The matrix $A$ is weakly diagonal dominant. If one of the rows of $W$ is trivial then it is even irreducibly diagonal dominant and by the Gershgorin circle theorem regular. Then $Ax=b$ has a solution.

If none of the rows of $W$ is trivial then the matrix is surely singular since the sum of the columns is zero. However the rank will be $n-1$, where $n$ is the size of the matrix. Indeed we may pick a column with at least one negative entry. If we remove this column and the corresponding row we obtain yet another irreducibly diagonal dominant matrix which is regular.

Hence the $Ax$ represents a hyperplane in $\mathbb R^n$ and the minimum of $J$ is the distance of $b$ to that hyperlane.

Edit: I am assuming that $A$ is irreducible. I actually have no idea whether this is ok in practise or whether restricting to "irreducible components" is a sensible thing to do.

Edit 2: If I am not mistaken there is (up to sign) exactly one normal vector $v$ orthogonal to the image of $A$ (that is in the $A$ singular case). Then the minimum of $J$ should be just the scalar product of $b$ and $v$. The vector $v$ should be orthogonal to all columns of $A$, which is the same as being in the kernel of $A^T$.

share|improve this answer
    
Good point. Actually, the trivial rows of $\mathbf W$ come from a constraint, fixing the value of certain $\mathbf x_i$. It makes sense that there is no solution without this constraint. Also, the property of diagonal dominance seems to be helpful to show that $\mathbf A$ is non-singular. I'll think about that! Thanks –  benni Jul 20 '12 at 13:56

Let's look at three cases of the problem $\mathbf{A}\mathbf{x} = \mathbf{b}$:

  1. $\mathbf{A}$ is square
  2. $\mathbf{A}$ is non-square and has row-dimension > column-dimension
  3. $\mathbf{A}$ is non-square and has column-dimension > row-dimension

Of course, these are the only three conditions that can happen.

WLOG, assume that the rows of $\mathbf{A}$ are linearly independent.

In case (1), the matrix is square and of full rank. Therefore, $\mathbf{A}$ is invertible, and you can exactly compute $\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$.

In case (2), you have an overdefined system. Think about it as if you've taken 4 measurements, with some possible measurement error, to solve for 2 unknowns. In this case, you cannot solve the system directly, but you can exactly compute the solution to $\mathbf{A}^T\mathbf{A}\mathbf{x} = \mathbf{A}^T\mathbf{b}$. This is exactly the same thing as computing a linear regression.

In case (3), you have an underdefined system. This is like trying to solve for 3 unknowns with only 2 equations. Very tricky indeed! In such a case, there is no direct solution to $\mathbf{A}^T\mathbf{A}\mathbf{x} = \mathbf{A}^T\mathbf{b}$. You have to, as in your problem, minimize the error using some optimization algorithm, of which there are many and varied.

How does this apply to your problem?

Well, you can completely determine your algorithmic approach by looking at the dimensions of the column space and row space of your matrix $\mathbf{A}$. By computing the number of linearly independent row/column vectors, you can determine easily whether $A$ represents a solvable, underdetermined, or overdetermined system, and choose your algorithmic approach based on that.

share|improve this answer
    
Thanks for your reply. Since, in my case, $\mathbf A$ is a square matrix, I obviously have to show that $\mathbf A$ is non-singular. Then $\mathbf A^T \mathbf A \mathbf x = \mathbf A^T \mathbf b$ should have the same unique solution as $\mathbf A \mathbf x = \mathbf b$. Right? –  benni Jul 20 '12 at 14:01
    
If $\mathbf{A}$ is non-singular, then $\mathbf{A}^T\mathbf{A}\mathbf{x} = \mathbf{A}^T\mathbf{b}$ will indeed have the same solution. However, it's a lot of computations, and for large $n$ that not only spends computing time unnecessarily, but may introduce error, leading to slightly different solutions. –  Arkamis Jul 20 '12 at 14:23
    
@benni To answer your final question in your OP, you can determine when $\mathbf{A}$ is invertible by checking to see if the rank of the matrix equals its square dimension, $n$. –  Arkamis Jul 20 '12 at 14:29
    
Yeah, right. However, I don't want to check this in practice, but show it in theory, so that I can solve $\mathbf A \mathbf x = \mathbf b$ instead of doing the costly matrix multiplication. –  benni Jul 20 '12 at 14:56
    
Well, if $\mathbf{A}$ is non-singular, then the solution $\mathbf{x}$ solving $\mathbf{Ax} = \mathbf{b}$ is unique. Then obviously $\mathbf{A}^T\mathbf{Ax} = \mathbf{A}^T\mathbf{b}$, because you're left-multiplying both sides of the equation by the same thing. Conversely, assume by contradiction that the solution to $\mathbf{A}^T\mathbf{Ax} = \mathbf{A}^T\mathbf{b}$ is $\mathbf{x'}$, different than the solution to $\mathbf{Ax} = \mathbf{b}$. Then, $\mathbf{A}^T\mathbf{Ax'} = \mathbf{A}^T\mathbf{b}$, and you left-multiply by $(\mathbf{A}^T)^{-1}$, obtaining $\mathbf{Ax'}=\mathbf{b}$. –  Arkamis Jul 20 '12 at 15:17

First two trivial remarks: any solution of $\mathbf A\mathbf x= \mathbf b$ is of course also a solution of $\mathbf A^T \mathbf A \mathbf x = \mathbf A^T\mathbf b $, and if $\mathbf A$ is invertible then the converse is also true (beacuse $\mathbf A^T$ is then also invertible and its inverse can be applied to both sides of the equation).

More is true though. Whether $\mathbf A$ is invertible or not, and even if it is non-square: if $\mathbf A\mathbf x= \mathbf b$ has any solutions, then $\mathbf A^T \mathbf A \mathbf x = \mathbf A^T\mathbf b $ has exactly the same solutions. However if $\mathbf A\mathbf x= \mathbf b$ has no solution then $\mathbf A^T \mathbf A \mathbf x = \mathbf A^T\mathbf b $ will still have solutions; a trivial example is when $\mathbf A$ is zero but $\mathbf b$ is not.

Here is the precise description of the situation. Let $c$ be the length of the columns of $\mathbf A$ (and of $\mathbf b$) and let $r\leq c$ be the rank of $\mathbf A$, which is the dimension of the image $I$ of the map $\mathbf x\mapsto\mathbf A\mathbf x$. By the rank-nullity theorem (and the fact that the rank of $\mathbf A^T$ is also $r$), the kernel $K$ of $\mathbf y\mapsto\mathbf A^t\mathbf y$ has co-dimension $r$ in $\mathbb R^c$ (and so dimension $c-r$). A crucial observation is that the subspaces $I,K$ of $\mathbb R^c$ form a direct sum (and $I\oplus K=\mathbb R^c$, due to the dimensions just mentioned): to be in $I$ means being a linear combination of the columns of $A$, and to be in $K$ means being orthogonal to all columns of $A$, and the two are compatible only for the zero vector. Therefore $(\mathbf b+K)\cap I$ consists of a single vector $\mathbf b'\in\mathbb R^c$, and the solutions of $\mathbf A^T \mathbf A \mathbf x = \mathbf A^T\mathbf b$ are those of $\mathbf A \mathbf x = \mathbf b'$.

Now if $\mathbf b=\mathbf b'$ this shows that the two equations of the question are equivalent, and this situation arises if and only if $\mathbf b\in I$. However if $\mathbf b\neq\mathbf b'$, so $\mathbf b\notin I$ and $\mathbf A \mathbf x = \mathbf b$ has no solution, then $\mathbf A \mathbf x = \mathbf b'$ still has solutions because $\mathbf b'\in I$ by construction. If $d$ is the length of the rows of $A$, the dimension of the set of solutions of $\mathbf A^T \mathbf A \mathbf x = \mathbf A^T\mathbf b$ is $d-r$ in all cases.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.