org.ejml.dense.row.mult.VectorVectorMult_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.mult;
import org.ejml.data.DMatrix1Row;
import org.ejml.data.DMatrixD1;
import org.ejml.data.DMatrixRMaj;
/**
* Operations that involve multiplication of two vectors.
*
* @author Peter Abeles
*/
public class VectorVectorMult_DDRM {
/**
*
* 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(DMatrixD1 x, DMatrixD1 y )
{
int m = x.getNumElements();
double total = 0;
for( int i = 0; i < m; i++ ) {
total += x.get(i) * y.get(i);
}
return total;
}
/**
*
* return = xT*A*y
*
*
* @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(DMatrixD1 x, DMatrixD1 A , DMatrixD1 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(DMatrixD1 x, DMatrixD1 A , DMatrixD1 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(DMatrixD1 x, DMatrixD1 y, DMatrix1Row 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 , DMatrixD1 x, DMatrixD1 y, DMatrix1Row 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,
DMatrixD1 u ,
DMatrixD1 x , DMatrixD1 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,
DMatrixRMaj A ,
DMatrixRMaj u , DMatrixRMaj w ,
DMatrixRMaj 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,
DMatrixRMaj A ,
DMatrixRMaj u ,
DMatrixRMaj 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];
}
}
}
}