org.ejml.alg.dense.mult.VectorVectorMult 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
* Computes the inner product of the two vectors. In geometry this is known as the dot product.
*
* ∑k=1:n xk * yk
* where x and y are vectors with n elements.
*
*
*
* These functions are often used inside of highly optimized code and therefor sanity checks are
* kept to a minimum. It is not recommended that any of these functions be used directly.
*
*
* @param x A vector with n elements. Not modified.
* @param y A vector with n elements. Not modified.
* @return The inner product of the two vectors.
*/
public static double innerProd( D1Matrix64F x, D1Matrix64F y )
{
int m = x.getNumElements();
double total = 0;
for( int i = 0; i < m; i++ ) {
total += x.get(i) * y.get(i);
}
return total;
}
/**
*
* xTAy
*
*
* @param x A vector with n elements. Not modified.
* @param A A matrix with n by m elements. Not modified.
* @param y A vector with m elements. Not modified.
* @return The results.
*/
public static double innerProdA( D1Matrix64F x, D1Matrix64F A , D1Matrix64F y )
{
int n = A.numRows;
int m = A.numCols;
if( x.getNumElements() != n )
throw new IllegalArgumentException("Unexpected number of elements in x");
if( y.getNumElements() != m )
throw new IllegalArgumentException("Unexpected number of elements in y");
double result = 0;
for( int i = 0; i < m; i++ ) {
double total = 0;
for( int j = 0; j < n; j++ ) {
total += x.get(j)*A.unsafe_get(j,i);
}
result += total*y.get(i);
}
return result;
}
/**
*
* xTATy
*
*
* @param x A vector with n elements. Not modified.
* @param A A matrix with n by n elements. Not modified.
* @param y A vector with n elements. Not modified.
* @return The results.
*/
// TODO better name for this
public static double innerProdTranA( D1Matrix64F x, D1Matrix64F A , D1Matrix64F y )
{
int n = A.numRows;
if( n != A.numCols)
throw new IllegalArgumentException("A must be square");
if( x.getNumElements() != n )
throw new IllegalArgumentException("Unexpected number of elements in x");
if( y.getNumElements() != n )
throw new IllegalArgumentException("Unexpected number of elements in y");
double result = 0;
for( int i = 0; i < n; i++ ) {
double total = 0;
for( int j = 0; j < n; j++ ) {
total += x.get(j)*A.unsafe_get(i,j);
}
result += total*y.get(i);
}
return result;
}
/**
*
* Sets A ∈ ℜ m × n equal to an outer product multiplication of the two
* vectors. This is also known as a rank-1 operation.
*
* A = x * y'
* where x ∈ ℜ m and y ∈ ℜ n are vectors.
*
*
* Which is equivalent to: Aij = xi*yj
*
*
*
* These functions are often used inside of highly optimized code and therefor sanity checks are
* kept to a minimum. It is not recommended that any of these functions be used directly.
*
*
* @param x A vector with m elements. Not modified.
* @param y A vector with n elements. Not modified.
* @param A A Matrix with m by n elements. Modified.
*/
public static void outerProd( D1Matrix64F x, D1Matrix64F y, RowD1Matrix64F A ) {
int m = A.numRows;
int n = A.numCols;
int index = 0;
for( int i = 0; i < m; i++ ) {
double xdat = x.get(i);
for( int j = 0; j < n; j++ ) {
A.set(index++ , xdat*y.get(j) );
}
}
}
/**
*
* Adds to A ∈ ℜ m × n the results of an outer product multiplication
* of the two vectors. This is also known as a rank 1 update.
*
* A = A + γ x * yT
* where x ∈ ℜ m and y ∈ ℜ n are vectors.
*
*
* Which is equivalent to: Aij = Aij + γ xi*yj
*
*
*
* These functions are often used inside of highly optimized code and therefor sanity checks are
* kept to a minimum. It is not recommended that any of these functions be used directly.
*
*
* @param gamma A multiplication factor for the outer product.
* @param x A vector with m elements. Not modified.
* @param y A vector with n elements. Not modified.
* @param A A Matrix with m by n elements. Modified.
*/
public static void addOuterProd( double gamma , D1Matrix64F x, D1Matrix64F y, RowD1Matrix64F A ) {
int m = A.numRows;
int n = A.numCols;
int index = 0;
if( gamma == 1.0 ) {
for( int i = 0; i < m; i++ ) {
double xdat = x.get(i);
for( int j = 0; j < n; j++ ) {
A.plus( index++ , xdat*y.get(j) );
}
}
} else {
for( int i = 0; i < m; i++ ) {
double xdat = x.get(i);
for( int j = 0; j < n; j++ ) {
A.plus( index++ , gamma*xdat*y.get(j));
}
}
}
}
/**
*
* Multiplies a householder reflection against a vector:
*
* y = (I + γ u uT)x
*
*
* The Householder reflection is used in some implementations of QR decomposition.
*
* @param u A vector. Not modified.
* @param x a vector. Not modified.
* @param y Vector where the result are written to.
*/
public static void householder( double gamma,
D1Matrix64F u ,
D1Matrix64F x , D1Matrix64F y )
{
int n = u.getNumElements();
double sum = 0;
for( int i = 0; i < n; i++ ) {
sum += u.get(i)*x.get(i);
}
for( int i = 0; i < n; i++ ) {
y.set( i , x.get(i) + gamma*u.get(i)*sum);
}
}
/**
*
* Performs a rank one update on matrix A using vectors u and w. The results are stored in B.
*
* B = A + γ u wT
*
*
* This is called a rank1 update because the matrix u wT has a rank of 1. Both A and B
* can be the same matrix instance, but there is a special rank1Update for that.
*
*
* @param gamma A scalar.
* @param A A m by m matrix. Not modified.
* @param u A vector with m elements. Not modified.
* @param w A vector with m elements. Not modified.
* @param B A m by m matrix where the results are stored. Modified.
*/
public static void rank1Update( double gamma,
DenseMatrix64F A ,
DenseMatrix64F u , DenseMatrix64F w ,
DenseMatrix64F B )
{
int n = u.getNumElements();
int matrixIndex = 0;
for( int i = 0; i < n; i++ ) {
double elementU = u.data[i];
for( int j = 0; j < n; j++ , matrixIndex++) {
B.data[matrixIndex] = A.data[matrixIndex] + gamma*elementU*w.data[j];
}
}
}
/**
*
* Performs a rank one update on matrix A using vectors u and w. The results are stored in A.
*
* A = A + γ u wT
*
*
* This is called a rank1 update because the matrix u wT has a rank of 1.
*
*
* @param gamma A scalar.
* @param A A m by m matrix. Modified.
* @param u A vector with m elements. Not modified.
*/
public static void rank1Update( double gamma,
DenseMatrix64F A ,
DenseMatrix64F u ,
DenseMatrix64F w )
{
int n = u.getNumElements();
int matrixIndex = 0;
for( int i = 0; i < n; i++ ) {
double elementU = u.data[i];
for( int j = 0; j < n; j++ ) {
A.data[matrixIndex++] += gamma*elementU*w.data[j];
}
}
}
}