org.ejml.ops.NormOps Maven / Gradle / Ivy
Show all versions of ejml Show documentation
/*
* Copyright (c) 2009-2012, Peter Abeles. All Rights Reserved.
*
* This file is part of Efficient Java Matrix Library (EJML).
*
* EJML is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation, either version 3
* of the License, or (at your option) any later version.
*
* EJML is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with EJML. If not, see
* Norms are a measure of the size of a vector or a matrix. One typical application is in error analysis.
*
*
* Vector norms have the following properties:
*
* - ||x|| > 0 if x ≠ 0 and ||0|| = 0
* - ||αx|| = |α| ||x||
* - ||x+y|| ≤ ||x|| + ||y||
*
*
*
*
* Matrix norms have the following properties:
*
* - ||A|| > 0 if A ≠ 0 where A ∈ ℜ m × n
* - || α A || = |α| ||A|| where A ∈ ℜ m × n
* - ||A+B|| ≤ ||A|| + ||B|| where A and B are ∈ ℜ m × n
* - ||AB|| ≤ ||A|| ||B|| where A and B are ∈ ℜ m × m
*
* Note that the last item in the list only applies to square matrices.
*
*
*
* Matrix norms can be induced from vector norms as is shown below:
*
* ||A||M = maxx≠0||Ax||v/||x||v
*
* where ||.||M is the induced matrix norm for the vector norm ||.||v.
*
*
*
* By default implementations that try to mitigate overflow/underflow are used. If the word fast is
* found before a function's name that means it does not mitigate those issues, but runs a bit faster.
*
*
* @author Peter Abeles
*/
public class NormOps {
/**
* Normalizes the matrix such that the Frobenius norm is equal to one.
*
* @param A The matrix that is to be normalized.
*/
public static void normalizeF( DenseMatrix64F A ) {
double val = normF(A);
if( val == 0 )
return;
int size = A.getNumElements();
for( int i = 0; i < size; i++) {
A.div(i , val);
}
}
/**
*
* The condition number of a matrix is used to measure the sensitivity of the linear
* system Ax=b. A value near one indicates that it is a well conditioned matrix.
*
* κp = ||A||p||A-1||p
*
*
* If the matrix is not square then the condition of either ATA or AAT is computed.
*
* @param A The matrix.
* @param p p-norm
* @return The condition number.
*/
public static double conditionP( DenseMatrix64F A , double p )
{
if( p == 2 ) {
return conditionP2(A);
} else if( A.numRows == A.numCols ){
// square matrices are the typical case
DenseMatrix64F A_inv = new DenseMatrix64F(A.numRows,A.numCols);
if( !CommonOps.invert(A,A_inv) )
throw new IllegalArgumentException("A can't be inverted.");
return normP(A,p) * normP(A_inv,p);
} else {
DenseMatrix64F pinv = new DenseMatrix64F(A.numCols,A.numRows);
CommonOps.pinv(A,pinv);
return normP(A,p) * normP(pinv,p);
}
}
/**
*
* The condition p = 2 number of a matrix is used to measure the sensitivity of the linear
* system Ax=b. A value near one indicates that it is a well conditioned matrix.
*
* κ2 = ||A||2||A-1||2
*
*
* This is also known as the spectral condition number.
*
*
* @param A The matrix.
* @return The condition number.
*/
public static double conditionP2( DenseMatrix64F A )
{
SingularValueDecomposition svd = DecompositionFactory.svd(A.numRows,A.numCols,false,false,true);
svd.decompose(A);
double[] singularValues = svd.getSingularValues();
int n = SingularOps.rank(svd,1e-12);
if( n == 0 ) return 0;
double smallest = Double.MAX_VALUE;
double largest = Double.MIN_VALUE;
for( double s : singularValues ) {
if( s < smallest )
smallest = s;
if( s > largest )
largest = s;
}
return largest/smallest;
}
/**
*
* This implementation of the Frobenius norm is a straight forward implementation and can
* be susceptible for overflow/underflow issues. A more resilient implementation is
* {@link #normF}.
*
*
* @param a The matrix whose norm is computed. Not modified.
*/
public static double fastNormF( D1Matrix64F a ) {
double total = 0;
int size = a.getNumElements();
for( int i = 0; i < size; i++ ) {
double val = a.get(i);
total += val*val;
}
return Math.sqrt(total);
}
/**
*
* Computes the Frobenius matrix norm:
*
* normF = Sqrt{ ∑i=1:m ∑j=1:n { aij2} }
*
*
* This is equivalent to the element wise p=2 norm. See {@link #fastNormF} for another implementation
* that is faster, but more prone to underflow/overflow errors.
*
*
* @param a The matrix whose norm is computed. Not modified.
* @return The norm's value.
*/
public static double normF( D1Matrix64F a ) {
double total = 0;
double scale = CommonOps.elementMaxAbs(a);
if( scale == 0.0 )
return 0.0;
final int size = a.getNumElements();
for( int i = 0; i < size; i++ ) {
double val = a.get(i)/scale;
total += val*val;
}
return scale*Math.sqrt(total);
}
/**
*
* Element wise p-norm:
*
* norm = {∑i=1:m ∑j=1:n { |aij|p}}1/p
*
*
*
* This is not the same as the induced p-norm used on matrices, but is the same as the vector p-norm.
*
*
* @param A Matrix. Not modified.
* @param p p value.
* @return The norm's value.
*/
public static double elementP( RowD1Matrix64F A , double p ) {
if( p == 1 ) {
return CommonOps.elementSumAbs(A);
} if( p == 2 ) {
return normF(A);
} else {
double max = CommonOps.elementMaxAbs(A);
if( max == 0.0 )
return 0.0;
double total = 0;
int size = A.getNumElements();
for( int i = 0; i < size; i++ ) {
double a = A.get(i)/max;
total += Math.pow(Math.abs(a),p);
}
return max*Math.pow(total,1.0/p);
}
}
/**
* Same as {@link #elementP} but runs faster by not mitigating overflow/underflow related problems.
*
* @param A Matrix. Not modified.
* @param p p value.
* @return The norm's value.
*/
public static double fastElementP( D1Matrix64F A , double p ) {
if( p == 2 ) {
return fastNormF(A);
} else {
double total = 0;
int size = A.getNumElements();
for( int i = 0; i < size; i++ ) {
double a = A.get(i);
total += Math.pow(Math.abs(a),p);
}
return Math.pow(total,1.0/p);
}
}
/**
* Computes either the vector p-norm or the induced matrix p-norm depending on A
* being a vector or a matrix respectively.
*
* @param A Vector or matrix whose norm is to be computed.
* @param p The p value of the p-norm.
* @return The computed norm.
*/
public static double normP( DenseMatrix64F A , double p ) {
if( p == 1 ) {
return normP1(A);
} else if( p == 2 ) {
return normP2(A);
} else if( Double.isInfinite(p)) {
return normPInf(A);
}
if( MatrixFeatures.isVector(A) ) {
return elementP(A,p);
} else {
throw new IllegalArgumentException("Doesn't support induced norms yet.");
}
}
/**
* An unsafe but faster version of {@link #normP} that calls routines which are faster
* but more prone to overflow/underflow problems.
*
* @param A Vector or matrix whose norm is to be computed.
* @param p The p value of the p-norm.
* @return The computed norm.
*/
public static double fastNormP( DenseMatrix64F A , double p ) {
if( p == 1 ) {
return normP1(A);
} else if( p == 2 ) {
return fastNormP2(A);
} else if( Double.isInfinite(p)) {
return normPInf(A);
}
if( MatrixFeatures.isVector(A) ) {
return fastElementP(A,p);
} else {
throw new IllegalArgumentException("Doesn't support induced norms yet.");
}
}
/**
* Computes the p=1 norm. If A is a matrix then the induced norm is computed.
*
* @param A Matrix or vector.
* @return The norm.
*/
public static double normP1( DenseMatrix64F A ) {
if( MatrixFeatures.isVector(A)) {
return CommonOps.elementSumAbs(A);
} else {
return inducedP1(A);
}
}
/**
* Computes the p=2 norm. If A is a matrix then the induced norm is computed.
*
* @param A Matrix or vector.
* @return The norm.
*/
public static double normP2( DenseMatrix64F A ) {
if( MatrixFeatures.isVector(A)) {
return normF(A);
} else {
return inducedP2(A);
}
}
/**
* Computes the p=2 norm. If A is a matrix then the induced norm is computed. This
* implementation is faster, but more prone to buffer overflow or underflow problems.
*
* @param A Matrix or vector.
* @return The norm.
*/
public static double fastNormP2( DenseMatrix64F A ) {
if( MatrixFeatures.isVector(A)) {
return fastNormF(A);
} else {
return inducedP2(A);
}
}
/**
* Computes the p=∞ norm. If A is a matrix then the induced norm is computed.
*
* @param A Matrix or vector.
* @return The norm.
*/
public static double normPInf( DenseMatrix64F A ) {
if( MatrixFeatures.isVector(A)) {
return CommonOps.elementMaxAbs(A);
} else {
return inducedPInf(A);
}
}
/**
*
* Computes the induced p = 1 matrix norm.
*
* ||A||1= max(j=1 to n; sum(i=1 to m; |aij|))
*
*
* @param A Matrix. Not modified.
* @return The norm.
*/
public static double inducedP1( DenseMatrix64F A ) {
double max = 0;
int m = A.numRows;
int n = A.numCols;
for( int j = 0; j < n; j++ ) {
double total = 0;
for( int i = 0; i < m; i++ ) {
total += Math.abs(A.get(i,j));
}
if( total > max ) {
max = total;
}
}
return max;
}
/**
*
* Computes the induced p = 2 matrix norm, which is the largest singular value.
*
*
* @param A Matrix. Not modified.
* @return The norm.
*/
public static double inducedP2( DenseMatrix64F A ) {
SingularValueDecomposition svd = DecompositionFactory.svd(A.numRows,A.numCols,false,false,true);
if( !svd.decompose(A) )
throw new RuntimeException("Decomposition failed");
double[] singularValues = svd.getSingularValues();
// the largest singular value is the induced p2 norm
return UtilEjml.max(singularValues,0,singularValues.length);
}
/**
*
* Induced matrix p = infinity norm.
*
* ||A||∞ = max(i=1 to m; sum(j=1 to n; |aij|))
*
*
* @param A A matrix.
* @return the norm.
*/
public static double inducedPInf( DenseMatrix64F A ) {
double max = 0;
int m = A.numRows;
int n = A.numCols;
for( int i = 0; i < m; i++ ) {
double total = 0;
for( int j = 0; j < n; j++ ) {
total += Math.abs(A.get(i,j));
}
if( total > max ) {
max = total;
}
}
return max;
}
}