org.ejml.dense.row.NormOps_DDRM Maven / Gradle / Ivy
/*
* Copyright (c) 2009-2017, Peter Abeles. All Rights Reserved.
*
* This file is part of Efficient Java Matrix Library (EJML).
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.ejml.dense.row;
import org.ejml.UtilEjml;
import org.ejml.data.DMatrix1Row;
import org.ejml.data.DMatrixD1;
import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.factory.DecompositionFactory_DDRM;
import org.ejml.interfaces.decomposition.SingularValueDecomposition_F64;
/**
*
* 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_DDRM {
/**
* 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( DMatrixRMaj 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(DMatrixRMaj A , double p )
{
if( p == 2 ) {
return conditionP2(A);
} else if( A.numRows == A.numCols ){
// square matrices are the typical case
DMatrixRMaj A_inv = new DMatrixRMaj(A.numRows,A.numCols);
if( !CommonOps_DDRM.invert(A,A_inv) )
throw new IllegalArgumentException("A can't be inverted.");
return normP(A,p) * normP(A_inv,p);
} else {
DMatrixRMaj pinv = new DMatrixRMaj(A.numCols,A.numRows);
CommonOps_DDRM.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( DMatrixRMaj A )
{
SingularValueDecomposition_F64 svd = DecompositionFactory_DDRM.svd(A.numRows,A.numCols,false,false,true);
svd.decompose(A);
double[] singularValues = svd.getSingularValues();
int n = SingularOps_DDRM.rank(svd,UtilEjml.TEST_F64);
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( DMatrixD1 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( DMatrixD1 a ) {
double total = 0;
double scale = CommonOps_DDRM.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(DMatrix1Row A , double p ) {
if( p == 1 ) {
return CommonOps_DDRM.elementSumAbs(A);
} if( p == 2 ) {
return normF(A);
} else {
double max = CommonOps_DDRM.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(DMatrixD1 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(DMatrixRMaj 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_DDRM.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(DMatrixRMaj 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_DDRM.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( DMatrixRMaj A ) {
if( MatrixFeatures_DDRM.isVector(A)) {
return CommonOps_DDRM.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( DMatrixRMaj A ) {
if( MatrixFeatures_DDRM.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( DMatrixRMaj A ) {
if( MatrixFeatures_DDRM.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( DMatrixRMaj A ) {
if( MatrixFeatures_DDRM.isVector(A)) {
return CommonOps_DDRM.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( DMatrixRMaj 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( DMatrixRMaj A ) {
SingularValueDecomposition_F64 svd = DecompositionFactory_DDRM.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( DMatrixRMaj 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;
}
}